This online documentation describes the necessary steps to perform the network optimization introduced in the paper What is the Minimal Systemic Risk in Financial Exposure Networks? by Diem, Pichler & Thurner (2019). A preprint of the paper is available on arXiv:1905.05931. The goal of this paper is to show the huge systemic risk reduction potential stemming from the network structure of financial exposure networks. We showed that for ten snap shots of the Austrian interbank market the proposed optimization substantially reduces systemic risk – measured with DebtRank (Battiston et al.,2012). First, we brefiefly outline the structure of the optimization problem and its formulation as mixed integer linear program (MILP). For the details of the reformulation and interpretation of the quantities we refer to the paper.

First, we construct a randomly generated sample network as reference network. Second, we provide the R code to create the objective function and the constraints to solve optimization problem with the R package ROI (Theußl et al., 2019). Finally we visualize the results for the randomly construced sample network.

# 1 The Network optimization Problem

Given an interbank network, $$L$$, we want to minimize its DebtRank. Because the calculation of DebtRank is an interative procedure it is difficult to optimize. Thus, we approximate DebtRank by the so called direct impacts of a bank on its neighbours (adjacent banks in the network). If bank $$i$$ defaults the direct loss of bank $$j$$ is $$\min \Big( \frac{L_{ij}}{e_j},1 \Big)$$. In the objective function bank $$j$$ has an economic weight corresponding to its share of all interbank market loans, $$\frac{a_j}{\bar{L}}$$. During the optimization the total amount of loans a bank has borrowed from (row sums of $$L$$) or lent to (column sums $$L$$) other banks shall stay the same. Additionally, the optimization shall leave the sum of loans weighted by their credit riskiness $$\kappa$$ unchanged. We showed that this problem can be formulated as \begin{align} \min_{L \in \{ M : \; M \in \mathbb{R}_+^{n \times n}, \; M_{ii}=0 \} } & \sum_{i=1}^{n}\sum_{j=1}^{n} \min \Big( \frac{L_{ij}}{e_j},1 \Big) \frac{a_j}{\bar{L}} \\ \text{s.t.} \qquad & l_i = \sum_{j=1}^{n} L_{ij} \quad \forall i \tag{1}\\ & a_i = \sum_{j=1}^{n} L_{ji} \quad \forall i \nonumber, \\ & r_i = \sum_{j=1}^{n} L_{ji}\kappa_j. \end{align} Details on the interpretation and definition of the quantities involved can be found in the paper.

## 1.1 Equivalent Mixed Integer Linear Program

We showed how the above optimization problem (Eq.1) can be rewritten into a mixed integer linear program (MILP), as

$\begin{eqnarray} \label{eq:optim_milp} \min_{ \{z \; \in \; \mathbb{R}_+^{4N^2} \} } & & c^\top z & \\ \text{s.t.} \quad & A_1z & \leq & 0 \tag{2} \\ & A_2z & = & a \\ & A_3z & = & l \\ & A_4z & = & r. \end{eqnarray}$

The first $$2N^2$$ entrities of $$z$$ correspond to the entrities of the interbank liability matrix $$L$$ by $$L_{11}=z_1+z_2$$, $$L_{21}=z_3+z_4$$, $$\dots$$,the second set of $$2N^2$$ variables are binary variables, where $$z_{2N^2+k}=1$$, if $$z_{k}>0$$ and zero otherwise. $$A_1$$ contains the constraints on the variables $$z$$ to ensure the equivalence of (2) and (1). $$A_2$$ and $$A_3$$ represent the constraints that the row and column sums need stay constant during the optimization and $$A_4$$ implents the credit risk constraint. For implementing other constraints with economic interpretations see Section 3.1.

## 1.2 Generate a Reference Network to Optimize

Before we can apply our optimization procedure we need an “empirically” observed financial network we can optimize. The optimization without the credit risk constraint can also be performed with the row and column sums of $$L$$ and the equity vector $$e$$. To illustrate the optimization, we use a randomly generated network with 25 banks. First, we sample the interbank asset vector $$a$$ and use a perturbed version as the interbank liability vector $$l$$. We set the equity vector $$e$$ proportional to the mean of $$a$$ and $$l$$. The sample is constructed such that we have 2 large banks, 3 medium sized banks and 20 smaller banks. Then we use the R package systemicrisk from Gandy et al. (2016) to generate a sample network. The exact construction of the sample network, $$L$$, is shwon in the code chunk below.


set.seed(50)

n_b <- 2; n_m <- 3; n_s <- 20
a <- c(runif(n_b,6000,10000), runif(n_m,2000,6000), runif(n_s,500,2000))
l <- a + c(rnorm(n_b,0,2000),rnorm(n_m,0,700),rnorm(n_s,0,150))

n <- length(a)
a[n+1] <- if(sum(l)-sum(a) >0){sum(l)-sum(a)}else{0}
l[n+1] <- if(sum(a)-sum(l) >0){sum(a)-sum(l)}else{0}
n <- length(a)
E <- (a+l)/5

# use the systemicrisk package
lambdaprior=Model.fitness.genlambdaparprior(ratescale=20))

thin <-100

res <- sample_HierarchicalModel(l=l, a=a, model=mod, nsamples=100, thin=thin, silent=TRUE)

# rounding the variables to two digits (monetary units)
L <- round(res$L[],2) a <- colSums(L) l <- rowSums(L) E <- round(E,2) kappa <- round(l/E,2) For the illustration we choose the credit risk proxy $$\kappa$$ to be interbank leverage ratio of the bank; i.e. kappa <- l/E. In practice different credit risk proxies would be desirable. # 2 Creating & Solving the MILP ## 2.1 Create Objective Functions and Constraints for MILP The quantities $$c, A_1, A_2, A_3, A_4$$ from (Eq.2) are generated with the code chunk below. The details of the construction for the objective function can be found in Diem et al. (2019), (Eq.19), and the details for the constraint matrices are explained in Appendix A. library(slam) N <- 4*n^2 # Construct the objective function vector c slopes <- rep(a/E, each=n) c <- c(as.vector(rbind(slopes,rep(0,n*n))), rep(0, 2*n*n)) # Construct the four constraint matrices A1, A2, A3 & A4 # A1 triplet_1 <- matrix(numeric(3*8*n^2), ncol=3) a_l <- rep(a, each = n) e_l <- rep(E, each = n) l_l <- rep(l, n) u_b <- pmax( rep(0,n^2), pmin(a_l, l_l) - e_l) for(i in 1:(n^2)){ #row col value triplet_1[(i*8-7):(i*8), ] <- rbind(c(4*i-3, 2*i-1, 1), c(4*i-3, 2*i-1+2*n^2, -e_l[i]), c(4*i-2, 2*i-1, -1), c(4*i-2, 2*i+2*n^2, e_l[i]), c(4*i-1, 2*i, 1), c(4*i-1, 2*i+2*n^2, -u_b[i]), c(4*i, 2*i-1+2*n^2, -1), c(4*i, 2*i+2*n^2, 1)) } A1 <- simple_triplet_matrix(triplet_1[,1], triplet_1[,2], triplet_1[,3]) # A2 triplet_2 <- matrix(numeric(3*(2*n^2)), ncol=3) for(i in 1:n){ for(j in 1:n){ v <- if(i != j){1}else{0} triplet_2[((2*i-2)*n+2*j-1):((2*i-2)*n+2*j), ] <- rbind( c(i, (2*j-2)*n +2*i-1, v), c(i, (2*j-2)*n +2*i, v)) } } A2 <- cbind( simple_triplet_matrix(triplet_2[,1], triplet_2[,2], triplet_2[,3]), simple_triplet_zero_matrix(n, 2*n^2)) # A3 triplet_3 <- matrix(numeric(3*(2*n^2)), ncol=3) for(i in 1:n){ for(j in 1:n){ v <- if(i != j){1}else{0} triplet_3[((2*i-2)*n+2*j-1):((2*i-2)*n+2*j), ] <- rbind( c(i, (2*i-2)*n +2*j-1, v), c(i, (2*i-2)*n +2*j, v)) } } A3 <- cbind( simple_triplet_matrix(triplet_3[,1], triplet_3[,2], triplet_3[,3]), simple_triplet_zero_matrix(n, 2*n^2)) # A4 triplet_4 <- matrix(numeric(3*(2*n^2)), ncol=3) for(i in 1:n){ for(j in 1:n){ v <- if(i != j){1*kappa[j]}else{0} triplet_4[((2*i-2)*n+2*j-1):((2*i-2)*n+2*j), ] <- rbind( c(i, (2*i-2)*n +2*j-1, v), c(i, (2*i-2)*n +2*j, v)) } } A4 <- cbind( simple_triplet_matrix(triplet_4[,1], triplet_4[,2], triplet_4[,3]), simple_triplet_zero_matrix(n, 2*n^2)) # Merge the constraints together A <- rbind(A1, A2, A3, A4) # Check the dimensions A #> A 2782x2704 simple triplet matrix. # Create upper bounds for the variables ub <- numeric(n^2) for(i in 1:n){ # col index for(j in 1:n){ # row index ub[((i-1)*n)+j] <- max(0, min( a[i],l[j])-E[i]) } } # Create a vector of zeros and ones to ensure that the diagonal elements L_ii=0 dia_constr <- matrix(rep(1,N/2), ncol = 2*n) for(j in 1:n){ # col dia_constr[j, (2*(j-1)+1):(2*(j-1)+2)] <- c(0,0) } dia_constr <- as.vector(t(dia_constr)) # Create a long vector of equity entries E_long <- rep(E, each = n) # Create a vector so that y_i1 is constraint by E_long and y_i2 by ub u_long <- as.vector(rbind(E_long, ub)) # combine the upper bounds with the zero upper bounds u <- u_long*dia_constr # Create b, t, d for ROI OP() b <- c(rep(0,N), l, a, crossprod(L,kappa)) t <- c(rep("C",2*n^2), rep("B", 2*n^2)) d <- c(rep("<=", N),rep("==", 3*n)) ## 2.2 Solving the optimization with ROI Now we can create the ROI object with the function OP() and solve it with ROI_solve(). Then, we transform the solution vector again into a liability matrix. For the optimization we use the cplex solver, which was working best in our experience. Smaller problems could also be solved with the open source solver SYMPHONY. library(ROI) LMP_min <- OP(objective = L_objective(c), constraints = L_constraint(L = A, dir = d, rhs = b), types = t, bounds = V_bound(ui = 1:(N/2), ub = u, nobj=N), maximum = FALSE) as.character(Sys.time()) #>  "2020-01-08 13:40:07" LMP_min_sol <- ROI_solve(LMP_min, "cplex") as.character(Sys.time()) #>  "2020-01-08 13:43:34" LMP_min_sol #> Optimal solution found. #> The objective value is: 6.377855e+04 LMP_min_sol$objval
#>  63778.55
LMP_min_sol$objval/sum(a) # direct impacts #>  1.263229 sol_vec <- LMP_min_sol$solution[1:(N/2)] # extract the non binary solution variables
L_vec   <- sol_vec[rep(c(TRUE,FALSE),N/4)] + sol_vec[rep(c(FALSE,TRUE),N/4)]

L_minDI <- matrix(L_vec, nrow=n)

# check if the constraints are fulfilled
c(sum(round(rowSums(L_minDI), 8) == round(rowSums(L), 8))==n,
sum(round(colSums(L_minDI), 8) == round(colSums(L), 8))==n)
#>  TRUE TRUE

Note that we can find a network with maximal direct impacts, simply by setting maximum = TRUE in the OP() function. This yields a network with a substantially higher DebtRank than the original one.

# maximisation problem
LMP_max <- OP(objective = L_objective(c),
constraints = L_constraint(L = A,
dir = d,
rhs = b),
types = t,
bounds = V_bound(ui = 1:(N/2), ub = u, nobj=N),
maximum = TRUE)
as.character(Sys.time())
#>  "2020-01-08 13:43:34"

LMP_max_sol <- ROI_solve(LMP_max, "cplex")

as.character(Sys.time())
#>  "2020-01-08 13:43:34"

LMP_max_sol
#> Optimal solution found.
#> The objective value is: 1.274026e+05
LMP_max_sol$objval #>  127402.6 LMP_max_sol$objval/sum(a)
#>  2.523397

sol_vec <- LMP_max_sol\$solution[1:(N/2)] # extract solutions of non binary objectives
L_vec   <- sol_vec[rep(c(TRUE,FALSE),N/4)] + sol_vec[rep(c(FALSE,TRUE),N/4)]

L_maxDI <- matrix(L_vec, nrow=n)

# check if the constraints are fulfilled
c(sum(round(rowSums(L_maxDI), 8) == round(rowSums(L), 8))==n,
sum(round(colSums(L_maxDI), 8) == round(colSums(L), 8))==n)
#>  TRUE TRUE

# 3 Comparing Original and Optimized Networks

We compare the direct impacts and DebtRanks of the three network types in Table 1. The optimized network show substantially smaller and higher DebtRank as the reference network.

Table 1, Comparison of Network DebtRanks
Minimized Empirical Maximized
DebtRank 4.48 9.69 17.91
Direct impacts 1.26 2.37 2.52

Figure 1 shows the systemic risk profile of the reference network and the minimized network. The optimization reduces the systemic risk for most banks. Figure 1, Systemic Risk Profile

## 3.1 Visualizing the Networks

The three networks are visualized in Figure 21. The node colors indicate the level of DebtRank, the node sizes correspond to bank equity and the link widths to the size of the loan.