In [34]:
import numpy as np
import networkx as nx
import matplotlib.pyplot as plt
def simulate_nested_split_rd(n, eta, rho, lam, zeta,
max_iter=100, tol=1e-8):
"""
Simulate the fixed‐point (A, y) for given parameters.
Args:
n : number of firms
eta : array‐like, length‐n vector of productivities η_i
rho : float, benefit parameter
lam : float, spill‐out parameter λ
zeta : float, link cost ζ
max_iter: maximum number of updates to A,y
tol : convergence tolerance on max |y_new - y_old|
Returns:
A : (n×n) adjacency matrix of the network
y : length‐n vector of R&D efforts
"""
η = np.asarray(eta).reshape(n)
# initialize A to no links
A = np.zeros((n,n), dtype=int)
# initial guess for y: just the stand‐alone η
y = η.copy()
for it in range(max_iter):
# --- 1) given A, solve for y: (I - B) y = η ---
B = np.zeros((n,n))
# off‐diagonal entries
B += rho * A
B -= lam * (np.ones((n,n)) - np.eye(n))
# zero out diagonal
np.fill_diagonal(B, 0.0)
# solve (I - B) y = η
M = np.eye(n) - B
y_new = np.linalg.solve(M, η)
# --- 2) update adjacency A by threshold rule ---
thresh = zeta / rho
# outer product y_i * y_j
Yprod = np.outer(y_new, y_new)
A_new = (Yprod > thresh).astype(int)
np.fill_diagonal(A_new, 0) # no self‐links
# check convergence
if np.max(np.abs(y_new - y)) < tol and np.array_equal(A_new, A):
break
y = y_new
A = A_new
return A, y
def plot_adjacency(A, ax=None):
"""Plot adjacency matrix as a heatmap."""
if ax is None:
fig, ax = plt.subplots()
im = ax.matshow(A, cmap='Blues')
ax.set_xlabel('j')
ax.set_ylabel('i')
ax.set_xticks(np.arange(-0.5, A.shape[1], 1), minor=True)
ax.set_yticks(np.arange(-0.5, A.shape[0], 1), minor=True)
ax.grid(which='minor', color='gray', linestyle='-', linewidth=0.5)
ax.tick_params(which='minor', bottom=False, left=False)
return ax
if __name__ == "__main__":
# === parameters from Figure 4 ===
n = 10
lam = 0.06
rho = 0.05 # in paper, rho = 0.02
eta = [1.00, 0.71, 0.58, 0.50, 0.45,
0.41, 0.38, 0.35, 0.33, 0.32]
zetas = [0.005, 0.0075, 0.015]
# Simulate for each ζ and plot
fig, axes = plt.subplots(1, len(zetas), figsize=(12,4))
for ax, ζ in zip(axes, zetas):
A, y = simulate_nested_split_rd(n, eta, rho, lam, ζ)
ax = plot_adjacency(A, ax=ax)
ax.set_title(fr'$\zeta={ζ}$')
plt.suptitle("Stepwise Adjacency Matrices of the Nested‐Split Graph")
plt.tight_layout(rect=[0,0,1,0.95])
plt.show()
