Does one Gibbs Step on a cycle

Description

Execute one Gibbs step on a cycle keeping row and column sums fixed

Usage

ERE_step_cycle(r, c, L, lambda, p, eps = 1e-10)

Arguments

Value

no return value

Examples

L=matrix(rexp(9),nrow=3)
lambda <- matrix(0.5,nrow=3,ncol=3)
p <- matrix(0.7, nrow=3,ncol=3)
ERE_step_cycle(r=c(0,1),c=c(1,2),L=L,lambda=lambda,p=p)
ERE_step_cycle(r=c(0,1,2),c=c(0,1,2),L=L,lambda=lambda,p=p)
ERE_step_cycle(r=c(0,1,2),c=c(2,1,0),L=L,lambda=lambda,p=p)

Gibbs sampling step of a matrix in the ERE model

Description

The sampling is conditional on row and column sums and uses k-cycle steps. Then dimensions of L, lambda and p must match.

Usage

GibbsSteps_kcycle(L, lambda, p, it = 1000L, eps = 1e-10, debug = 0L)

Arguments

Value

no return value

Examples

L <- matrix(c(1,2,3,4,5,6,7,8,9),nrow=3)
diag(L) <- 0
lambda <- matrix(0.5,nrow=3,ncol=3)
p <- matrix(0.7, nrow=3,ncol=3)
diag(p) <- 0
GibbsSteps_kcycle(L=L,lambda=lambda,p=p)
L
L <- matrix(1:16,nrow=4)
diag(L) <- 0
lambda <- matrix(0.5,nrow=4,ncol=4)
p <- matrix(0.25, nrow=4,ncol=4)
diag(p) <- 0
GibbsSteps_kcycle(L=L,lambda=lambda,p=p)
L

Combination of Independent Models for p and lambda

Description

Combination of Independent Models for p and lambda

Usage

Model.Indep.p.lambda(model.p, model.lambda)

Arguments

Value

the resulting model.

Examples

n <- 5
m <- Model.Indep.p.lambda(Model.p.BetaPrior(n),
                          Model.lambda.GammaPrior(n,scale=1e-1))
genL(m)

Fitness model for liabilities matrix

Description

Assumes a diagonal consisting of 0s.

Usage

Model.additivelink.exponential.fitness(
  n,
  alpha,
  beta,
  gamma = 1,
  lambdaprior,
  sdpropfitness = 1/sqrt(n)
)

Arguments

Value

A model to be used by sample_HierarchicalModel. This is a list of functions. It includes a function accrates() that repors acceptance rates for the Metropolis-Hasting steps involved.

Examples

mod <- Model.additivelink.exponential.fitness(10,alpha=-2.5,beta=0.1,
                lambdaprior=Model.fitness.genlambdaparprior(ratescale=1e3))
theta <- mod$rtheta()
L <- genL(mod)
l <- rowSums(L$L)
a <- colSums(L$L)
## increase number of samples and thinning in real examples
res <- sample_HierarchicalModel(l=l,a=a,model=mod,nsamples=4,thin=50)
mod$accrates()

Mean out-degree of a node with given fitness in the fitness model

Description

Computes the mean out-degree of a node with given fitness x in the fitness model implemented in Model.additivelink.exponential.fitness. The function returns the mean out-degree divided by n-1.

Usage

Model.fitness.conditionalmeandegree(x, alpha, beta, gamma = 1)

Arguments

Value

Mean out-degree divided by n-1.

Prior distribution for eta and zeta in the fitness model

Description

Assumes a uniform distribution on the shape parameter zeta and an exponential distribution on the scale parameter eta. To be used as prior for Model.additivelink.exponential.fitness.

Usage

Model.fitness.genlambdaparprior(
  shapemin = 0.75,
  shapemax = 1.5,
  ratescale,
  sdshapeprob = 0.1,
  sdpropscale = 0.1
)

Arguments

Value

list of functions necessary for constructing Metropolis-Hastings updates.

Mean out-degree of a random node the fitness model

Description

Computes the relative mean out-degree of a randomly chosen node given fitness x in the fitness model implemented in Model.additivelink.exponential.fitness. The function returns the mean out-degree divided by n-1.

Usage

Model.fitness.meandegree(alpha, beta, gamma = 1)

Arguments

Value

Mean out-degree divided by n-1.

Model with Gamma Prior on Lambda

Description

Assumes that all elements of lambda are equal to a parameter \theta, which has a Gamma prior.

Usage

Model.lambda.GammaPrior(n, shape = 1, scale = 1)

Arguments

Value

the resulting model.

Model Using Multiple Independent Components

Description

Assumes a multivariate hyperparameter \theta with each component following an independent Beta distribution. A matrix indicates which component \theta is used for what component of lambda.

Usage

Model.lambda.Gammaprior_mult(Ilambda, shape = 1, scale = 1)

Arguments

Value

the resulting model.

Model for a Constant lambda

Description

This model assumes that the parameter lambda is known.

Usage

Model.lambda.constant(lambda, n)

Arguments

Value

the resulting model.

Examples

m <- Model.lambda.constant(n=5,lambda=0.25)
m$matr(m$rtheta())
lambda<-matrix(c(NA,1,1,1e-4,NA,1e-4,1e4,1e4,NA),nrow=3)
m <- Model.lambda.constant(n=3,lambda=lambda)
m$matr(m$rtheta())

Model for a Constant lambda and Non-Square Matrices

Description

This model assumes that the parameter lambda is known.

Usage

Model.lambda.constant.nonsquare(lambda, nrow, ncol)

Arguments

Value

the resulting model.

Examples

m <- Model.lambda.constant.nonsquare(nrow=5,ncol=3,lambda=0.25)
m$matr(m$rtheta())

Model for a Random One-dimensional p

Description

Assumes a Beta prior on the one-dimensional link existence probabilities p. This model has a one-dimensional parameter.

Usage

Model.p.BetaPrior(n, shape1 = 1, shape2 = 1)

Arguments

Value

the resulting model.

Examples

m <- Model.p.BetaPrior(5)
m$matr(m$rtheta())

Model Using Multiple Independent Components

Description

Assumes a multivariate hyperparameter \theta with each component following an independent Beta distribution. A matrix indicates which component \theta is used for what component of p.

Usage

Model.p.Betaprior_mult(Ip, shape1 = 1, shape2 = 1)

Arguments

Value

the resulting model.

Multiplicative Fitness Model for Power Law

Description

This model has a power law of the degree distribution with a parameter \alpha and is tuned to a desired link existence probability. It is based on a fitness model.

Usage

Model.p.Fitness.Servedio(n, alpha, meandegree, sdprop = 0.1)

Arguments

Details

Every node i has a fitness \theta_i being an independent realisation of a U[0,1] distribution. The probability of a link between a node with fitness x and a node with fitness y is g(x)g(y) where g is as follows. If \alpha=-1 then

g(x)=g0*\exp(-\log(g0)*x)

Otherwise,

g(x)=(g0^(\alpha+1)+(1-g0^(\alpha+1))*x)^(1/(\alpha+1))

where g0 is tuned numerically to achieve the desired overall mean degree.

Updating of the model parameters in the MCMC setup is done via a Metropolis-Hastings step, adding independent centered normal random variables to each node fitness in \theta.

Value

the resulting model.

References

Servedio V. D. P. and Caldarelli G. and Butta P. (2004) Vertex intrinsic fitness: How to produce arbitrary scale-free networks. Physical Review E 70, 056126.

Examples

n <- 5
mf <- Model.p.Fitness.Servedio(n=n,alpha=-2.5,meandegree=0.5)
m <- Model.Indep.p.lambda(model.p=mf,
                          model.lambda=Model.lambda.GammaPrior(n,scale=1e-1))
x <- genL(m)
l <- rowSums(x$L)
a <- colSums(x$L)
res <- sample_HierarchicalModel(l,a,model=m,nsamples=10,thin=10)

Model for a Constant p

Description

This model assumes that the link existence probabilities of the matrix are known.

Usage

Model.p.constant(n, p)

Arguments

Value

the resulting model.

Examples

m <- Model.p.constant(5,0.25)
m$matr(m$rtheta())

p <- matrix(c(0,0.99,0.99,0.5,0.5,0,0.01,0.01,0),nrow=3)
m <- Model.p.constant(5,p)
m$matr(m$rtheta())

Model for a constant p and Non-Square Matrices

Description

This model assumes that the link existence probabilities of the matrix are known.

Usage

Model.p.constant.nonsquare(nrow, ncol, p)

Arguments

Value

the resulting model.

Examples

m <- Model.p.constant.nonsquare(5,3,0.25)
m$matr(m$rtheta())

Calibrate ER model to a given density

Description

The model is an Erdos-Renyi model for the existence of links (a link exists independently of other links with a fixed probability) and the weight of each existing link follows an exponential distribution with a fixed rate parameter. This function chooses the two parameters such that the density of the network (the average proportion of existing links) is a certain desired value. Diagonal elements are being set to 0.

Usage

calibrate_ER(
  l,
  a,
  targetdensity,
  L_fixed = NA,
  nsamples_calib = 100,
  thin_calib = 100
)

Arguments

Value

Model that can be used to generate the desired samples using sample_HierarchicalModel.

Examples

## first generate a true network
n <- 10 # size of network
p <- 0.45
lambda <- 0.1
L <- matrix(nrow=n,rbinom(n*n,prob=p,size=1)*rexp(n*n,rate=lambda))

# then reconstruct with a target density of 0.55
model <- calibrate_ER(l=rowSums(L),a=colSums(L),
                      targetdensity=0.55,nsamples_calib=10)
Lsamp <- sample_HierarchicalModel(l=rowSums(L),a=colSums(L),model=model,
                                    nsamples=10,thin=1e2)
# check row sums
rowSums(L)
rowSums(Lsamp$L[[10]])
# check calibration
mean(Lsamp$L[[10]]>0)

# now an example with some fixed entries
L_fixed <- L
L_fixed[1:(n/2),] <- NA
# then reconstruct with a target density of 0.9
model <- calibrate_ER(l=rowSums(L),a=colSums(L),L_fixed=L_fixed,
                              targetdensity=0.9,nsamples_calib=10)
Lsamp <- sample_HierarchicalModel(l=rowSums(L),a=colSums(L),L_fixed=L_fixed,
                                  model=model,nsamples=10,thin=1e2)
mean(Lsamp$L[[10]][-(1:(n/2)),]>0) # known entries
mean(Lsamp$L[[10]][(1:(n/2)),]>0)  #reconstructed entries

Calibrate ER model to a given density with a nonsquare matrix

Description

The model is an Erdos-Renyi model for the existence of links (a link exists independently of other links with a fixed probability) and the weight of each existing link follows an exponential distribution with a fixed rate parameter. This function chooses the two parameters such that the density of the network (the average proportion of existing links) is a certain desired value. This function does not set diagonal values to 0.

Usage

calibrate_ER.nonsquare(
  l,
  a,
  targetdensity,
  L_fixed = NA,
  nsamples_calib = 100,
  thin_calib = 100
)

Arguments

Value

Model that can be used to generate the desired samples using sample_HierarchicalModel.

Examples

## first generate a true network
nrow <- 10 # size of network
ncol <- 8 # size of network
p <- 0.45
lambda <- 0.1
L <- matrix(nrow=nrow,rbinom(nrow*ncol,prob=p,size=1)*rexp(nrow*ncol,rate=lambda))

# then reconstruct with a target density of 0.55
model <- calibrate_ER.nonsquare(l=rowSums(L),a=colSums(L),
                      targetdensity=0.55,nsamples_calib=10)
Lsamp <- sample_HierarchicalModel(l=rowSums(L),a=colSums(L),model=model,
                                    nsamples=10,thin=1e2)
# check row sums
rowSums(L)
rowSums(Lsamp$L[[10]])
# check calibration
mean(Lsamp$L[[10]]>0)

# now an example with some fixed entries
L_fixed <- L
L_fixed[1:(nrow/2),] <- NA
# then reconstruct with a target density of 0.9
model <- calibrate_ER.nonsquare(l=rowSums(L),a=colSums(L),L_fixed=L_fixed,
                              targetdensity=0.9,nsamples_calib=10)
Lsamp <- sample_HierarchicalModel(l=rowSums(L),a=colSums(L),L_fixed=L_fixed,
                                  model=model,nsamples=10,thin=1e2)
mean(Lsamp$L[[10]][-(1:(nrow/2)),]>0) # known entries
mean(Lsamp$L[[10]][(1:(nrow/2)),]>0)  #reconstructed entries

Calibrate empirical fitness model to a given density

Description

The model is an empirical fitness based model for the existence of links (more details below) which contains one fixed parameter and the weight of each existing link follows an exponential distribution with a fixed rate parameter. This function chooses the two parameters such that the density of the network (the average proportion of existing links) with these given row and column sums is a certain desired value.

Usage

calibrate_FitnessEmp(
  l,
  a,
  targetdensity,
  L_fixed = NA,
  nsamples_calib = 100,
  thin_calib = 100
)

Arguments

Details

The empirical fitness model assumes that every node f[i]=log(l[i]+a[i]) has a fitness given by the observered row and column sum and that the existence probability of a link between node i and j is then given by p[i,j]=1/(1+exp(-(alpha+f[i]+f[j]))), where alpha is an additional parameter. The resulting model uses observed quantities (the row and column sums of the matrix) as input to the model and is thus an empirical Bayes approach.

Value

Model that can be used to generate the desired samples using sample_HierarchicalModel.

Examples

## first generate a true network
n <- 10 # size of network
ftrue <- rnorm(n) # vector of underlying fitnesses
p <- outer(ftrue,ftrue,FUN=function(x,y) 1/(1+exp(-(x+y))))
lambda <- 0.1
L <- matrix(nrow=n,rbinom(n*n,prob=p,size=1)*rexp(n*n,rate=lambda))

# then reconstruct with a target density of 0.7
model <- calibrate_FitnessEmp(l=rowSums(L),a=colSums(L),
                              targetdensity=0.7,nsamples_calib=10,thin_calib=50)
Lsamp <- sample_HierarchicalModel(l=rowSums(L),a=colSums(L),model=model,
                                    nsamples=10,thin=1e2)
# check row sums
rowSums(L)
rowSums(Lsamp$L[[10]])
# check calibration
mean(Lsamp$L[[10]]>0)

# now an example with some fixed entries
L_fixed <- L
L_fixed[1:(n/2),] <- NA
# then reconstruct with a target density of 0.9
model <- calibrate_FitnessEmp(l=rowSums(L),a=colSums(L),L_fixed=L_fixed,
                              targetdensity=0.9,nsamples_calib=10,thin_calib=50)
Lsamp <- sample_HierarchicalModel(l=rowSums(L),a=colSums(L),L_fixed=L_fixed,
                                  model=model,nsamples=10,thin=1e2)
mean(Lsamp$L[[10]][-(1:(n/2)),]>0) # known entries
mean(Lsamp$L[[10]][(1:(n/2)),]>0)  #reconstructed entries

Calibrate Thinning

Description

Attempts to automatically choose a thinning paramter to achieve an overall relative effective sample size (defined as the effective sample size divided by the number of samples) for all parameters in the model (that do not seem to be constant). This function provides no guarantees that the desired relative effective sample size (rESS) will actually be achieved - it is best treated as a rough guide for this.

Usage

choosethin(
  l,
  a,
  L_fixed = NA,
  model,
  relESStarget = 0.3,
  burnin = 100,
  matrpertheta = length(l)^2,
  silent = TRUE,
  maxthin = 10000
)

Arguments

Details

The approach used involves a pilot run of the sampler, followed by a computation of the acf (autocorrelation function) for each component. The acf is used only up to (and excluding) the point used where it becomes negative for the first time. This part of the acf is then used to approximate the rESS and to determine the amount of thinning needed. The reported result is the thinning needed to achieve the rESS for all components (the matrix as well as the parameter theta). The initial pilot run may not be sufficient and further pilot runs may have to be started.

Value

An integer describing the amount of thinning required.

Examples

set.seed(12689)
n <- 10
m <- Model.Indep.p.lambda(Model.p.BetaPrior(n),
                          Model.lambda.GammaPrior(n,scale=1e-1))
x <- genL(m)
l <- rowSums(x$L)
a <- colSums(x$L)
choosethin(l,a,model=m)
choosethin(l,a,model=m,relESStarget=0.7)

Creates a deep copy of a matrix

Description

Useful when calling ERE_step_cycle or GibbsSteps_kcycle to ensure that there are no side effects for the return values.

Usage

cloneMatrix(M)

Arguments

Value

A deep copy of the matrix.

Examples

lambda <- matrix(0.5,nrow=2,ncol=2)
p <- matrix(0.7, nrow=2,ncol=2)
L <- matrix(rexp(4),nrow=2);
L
Lold <- L
Lcopy <- cloneMatrix(L)
ERE_step_cycle(r=c(0,1),c=c(0,1),L=L,lambda=lambda,p=p)

L     ## new value
Lold  ## equal to L !!!
Lcopy ## still has the original value

Default of Banks

Description

Computes bank defaults based on a liabilities matrix and external assets and liabilities.

Usage

default(L, ea, el = 0, method = c("clearing", "cascade"), ...)

Arguments

Value

A list with at least one element "defaultind", which is a vector indicating which banks default (1=default, 0= no default). Depending on the method, other results such as the clearing vector may also be reported.

See Also

default_cascade, default_clearing,

Examples

ea <- c(1/2,5/8,3/4)
el <- c(3/2,1/2,1/2)
x <- 0.5
L <- matrix(c(0,x,1-x,1-x,0,x,x,1-x,0),nrow=3)
default(L,ea,el)
default(L,ea,el,"cascade")

Default Cascade

Description

Computes bank defaults via the default cascade algorithm.

Usage

default_cascade(L, ea, el = 0, recoveryrate = 0)

Arguments

Value

vector indicating which banks default (1=default, 0= no default)

Examples

ea <- c(1/2,5/8,3/4)
el <- c(3/2,1/2,1/2)
x <- 0.5
L <- matrix(c(0,x,1-x,1-x,0,x,x,1-x,0),nrow=3)
default_cascade(L,ea,el)

Clearing Vector with Bankruptcy Costs

Description

Computes bank defaults for the clearing vector approach without and with bankruptcy costs (Eisenberg and Noe, 2001), (Rogers and Veraart, 2013).

Usage

default_clearing(L, ea, el = 0, alpha = 1, beta = 1)

Arguments

Details

Without bankruptcy costs the approach of Eisenberg and Noe (2001) is used using a linear programme. With bankruptcy costs, the implementation is based on the Greatest Clearing Vector Algorithm (GA), see Definition 3.6, Rogers & Veraart (2013).

Value

A list consisting of a vector indicating which banks default (1=default, 0= no default) and the greatest clearing vector.

References

Eisenberg, L. and Noe, T.H. (2001). Systemic risk in financial systems. Management Science 47, 236–249.

Rogers, L. C. G. and Veraart, L. A. M. (2013) Failure and Rescue in an Interbank Network, Management Science 59 (4), 882–898.

Examples

ea <- c(1/2,5/8,3/4)
el <- c(3/2,1/2,1/2)
x <- 0.5
L <- matrix(c(0,x,1-x,1-x,0,x,x,1-x,0),nrow=3)
default_clearing(L,ea,el)
default_clearing(L,ea,el, alpha=0.5, beta=0.7)

Outputs Effective Sample Size Diagonistics for MCMC run

Description

Computes the Effective Sample Size using the method effectiveSize in of the package coda.

Usage

diagnose(res)

Arguments

Details

Currently only works with L where the diagonal is 0. The function ignores the diagonal and tries to determine from the row and column sums which parts of the matrix are 0.

Value

No return value. Called for printing the diagnostics.

Finds a Nonnegative Matrix Satisfying Row and Column Sums

Description

Given row and column sums and a matrix p which indicates which elements of the matrix can be present, this function computes a nonnegative matrix that match these row and column sums. If this is not possible then the function returns an error message.

Usage

findFeasibleMatrix(r, c, p, eps = 1e-09)

Arguments

Details

The function transforms the problem into a Maximum Flow problem of a graph and uses the Edmonds-Karps algorithm to solve it. If the error message "Could not find feasible matrix." is produced then this could be due to p imposing disconnected components in the graph implied by row and column sums that are not compatible with the row and column sums..

Value

A feasible matrix.

Examples

p=matrix(c(1,0,0,1),nrow=2)
findFeasibleMatrix(c(1,1),c(1,1),p=p)

n <- 4
M <- matrix(nrow=n,ncol=n,rexp(n*n)*(runif(n*n)>0.6))
M
r <- rowSums(M)
c <- colSums(M)
Mnew <- findFeasibleMatrix(r=r,c=c,p=(M>0)*0.5)
Mnew
rowSums(M);rowSums(Mnew)
colSums(M);colSums(Mnew)

Creates a feasible starting matrix with a desired mean average degree

Description

This extension of findFeasibleMatrix attempts to create a feasible matrix where a certain proportion of the entries is positive. There is no guarantee that this proportion is achieved. If it is not possible then this matrix will report a warning and simply return the matrix constructed by findFeasibleMatrix.

Usage

findFeasibleMatrix_targetmean(r, c, p, eps = 1e-09, targetmean = 0.3)

Arguments

Value

A feasible matrix.

Generate Liabilities Matrix from Prior

Description

Generates a libabilities matrix using a the prior distribution from a given model for p and lambda.

Usage

genL(model)

Arguments

Value

A list consisting of a liabilities matrix and the parameter vector theta.

Examples

n <- 5
m <- Model.Indep.p.lambda(Model.p.BetaPrior(n),
                          Model.lambda.GammaPrior(n,scale=1e-1))
genL(m)

Creates a feasible starting matrix

Description

Creates a matrix with nonnegative entries, given row and column sums and 0 on the diagonal. Superseeded by the more flexible findFeasibleMatrix.

Usage

getfeasibleMatr(L, A)

Arguments

Value

A matrix with nonnegative entries and given row/column sums and 0 on the diagonal.

Examples

getfeasibleMatr(c(0.5,1,0),c(0.5,0,1))
getfeasibleMatr(rep(1,4),rep(1,4))
getfeasibleMatr(2^(1:3),2^(3:1))
getfeasibleMatr(1:5,1:5)
getfeasibleMatr(1:5,5:1)

Sample from the ERE model with given row and column sums

Description

Samples from the Erdos Reny model with Exponential weights and known marginals. Runs a Gibbs sampler to do this. A starting liabilities is generated via getfeasibleMatr before steps_ERE is called.

Usage

sample_ERE(l, a, p, lambda, nsamples = 10000, thin = 1000, burnin = 10000)

Arguments

Value

List of simulation results

Examples

l <- c(1,2.5,3)
a <- c(0.7,2.7,3.1)
L <- sample_ERE(l,a,p=0.5,lambda=0.25,nsamples=5,thin=20,burnin=10)
L

Sample from Hierarchical Model with given Row and Column Sums

Description

Sample from Hierarchical Model with given Row and Column Sums

Usage

sample_HierarchicalModel(
  l,
  a,
  L_fixed = NA,
  model,
  nsamples = 10000,
  thin = choosethin(l = l, a = a, L_fixed = L_fixed, model = model, matrpertheta =
    matrpertheta, silent = silent),
  burnin = NA,
  matrpertheta = length(l)^2,
  silent = FALSE,
  tol = .Machine$double.eps^0.25
)

Arguments

Value

The resulting samples. A list with the first element, L, giving the samples of matrices, and the second element, theta, giving the samples of the hyperparameter (if hyperparameters are present).

Examples

n <- 10
m <- Model.Indep.p.lambda(Model.p.BetaPrior(n),
                          Model.lambda.GammaPrior(n,scale=1e-1))
x <- genL(m)
l <- rowSums(x$L)
a <- colSums(x$L)

res <- sample_HierarchicalModel(l,a,model=m)

# fixing one values
L_fixed <- matrix(NA,ncol=n,nrow=n)
L_fixed[1,2:5] <- x$L[1,2:5]

res <- sample_HierarchicalModel(l,a,model=m,L_fixed=L_fixed,
                                nsamples=1e2)
sapply(res$L,function(x)x[1,2:5])

Perform Steps of the Gibbs Sampler of the ERE model

Description

Runs a Gibbs sampler in the Erdos Reny model with Exponential weights (ERE model) and fixed marginals. The algorithm starts from a given matrix.

Usage

steps_ERE(L, p, lambda, nsamples = 10000, thin = 1000, burnin = 10000)

Arguments

Value

List of simulation results

See Also

sample_ERE

Examples

L <- matrix(rexp(4*4),nrow=4,ncol=4); diag(L)=0;
p <- matrix(0.5,nrow=4,ncol=4); diag(p) <-0;
lambda <- matrix(1,nrow=4,ncol=4); diag(lambda)<-0;

L <- steps_ERE(L=L,p=p,lambda=lambda,nsamples=5,thin=50,burnin=20)
L