The location quotient (LQ) measures industry $i$'s relative concentration in region $r$ compared to the national average:
$$LQ_{ir} = \frac{e_{ir}/E_r}{e_{in}/E_n}$$
where $e_{ir}$ is employment in industry $i$, region $r$; $E_r$ is total regional employment; $e_{in}$ and $E_n$ are corresponding national totals. $LQ > 1$ implies specialization; $LQ > 1.25$ is often used as an export-base threshold.
Shift-Share Analysis
The classic Dunn (1960) decomposition splits regional employment change $\Delta e_{ir}$ into three components:
import numpy as np
import pandas as pd
import scipy.linalg as la
import statsmodels.api as sm
import matplotlib.pyplot as plt
print("Regional economics environment ready")
Core Workflow
Step 1: Location Quotients and Shift-Share Analysis
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# -----------------------------------------------------------------# Build a stylized 4-sector IO table# Sectors: Agriculture, Manufacturing, Services, Construction# -----------------------------------------------------------------
sectors = ["Agriculture", "Manufacturing", "Services", "Construction"]
n = len(sectors)
# Intermediate transactions matrix Z (4x4)
Z = np.array([
[20, 80, 10, 5], # Agriculture sells to ...
[15, 120, 60, 30], # Manufacturing sells to ...
[ 5, 40, 150, 20], # Services sells to ...
[ 2, 15, 10, 8], # Construction sells to ...
], dtype=float)
# Final demand vector f
f = np.array([100, 200, 300, 80], dtype=float)
# Gross output X = Z.sum(axis=1) + f
X = Z.sum(axis=1) + f
print("Gross Output X:", X)
# -----------------------------------------------------------------# Technical coefficient matrix A# -----------------------------------------------------------------
A = Z / X[np.newaxis, :] # a_ij = z_ij / X_jprint("\nTechnical Coefficient Matrix A:")
print(pd.DataFrame(A, index=sectors, columns=sectors).round(3))
# -----------------------------------------------------------------# Leontief Inverse L = (I - A)^{-1}# -----------------------------------------------------------------
I = np.eye(n)
L = np.linalg.inv(I - A)
print("\nLeontief Inverse L:")
print(pd.DataFrame(L, index=sectors, columns=sectors).round(3))
# Output multipliers: column sums of L
output_mult = L.sum(axis=0)
print("\nOutput Multipliers (column sums of L):")
for s, m inzip(sectors, output_mult):
print(f" {s}: {m:.3f}")
# -----------------------------------------------------------------# Employment multipliers (if we have employment coefficients)# -----------------------------------------------------------------# Suppose employment per unit output (jobs per $M output)
emp_coeff = np.array([0.05, 0.03, 0.04, 0.06]) # l vector
emp_mult = emp_coeff @ L # total employment requirementsprint("\nEmployment Multipliers:")
for s, m inzip(sectors, emp_mult):
print(f" {s}: {m:.4f} jobs per $M final demand")
# -----------------------------------------------------------------# Simulate a demand shock: +$50M final demand in Manufacturing# -----------------------------------------------------------------
delta_f = np.array([0, 50, 0, 0], dtype=float)
delta_X = L @ delta_f
print(f"\nImpact of +$50M demand shock in Manufacturing:")
for s, dx inzip(sectors, delta_X):
print(f" {s}: +${dx:.1f}M output")
# -----------------------------------------------------------------# Visualization: Leontief inverse heatmap# -----------------------------------------------------------------
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
im = axes[0].imshow(L, cmap="Blues", aspect="auto")
axes[0].set_xticks(range(n)); axes[0].set_xticklabels(sectors, rotation=30, ha="right")
axes[0].set_yticks(range(n)); axes[0].set_yticklabels(sectors)
axes[0].set_title("Leontief Inverse (I-A)⁻¹")
plt.colorbar(im, ax=axes[0])
# Annotate cellsfor i inrange(n):
for j inrange(n):
axes[0].text(j, i, f"{L[i,j]:.2f}", ha="center", va="center",
color="white"if L[i,j] > 1.5else"black", fontsize=8)
# Output and employment multipliers bar chart
x = np.arange(n)
w = 0.35
axes[1].bar(x - w/2, output_mult, w, label="Output Multiplier", color="steelblue")
axes[1].bar(x + w/2, emp_mult / emp_mult.max() * output_mult.max(), w,
label="Emp. Multiplier (scaled)", color="orange")
axes[1].set_xticks(x); axes[1].set_xticklabels(sectors, rotation=20, ha="right")
axes[1].set_title("Output and Employment Multipliers")
axes[1].set_ylabel("Multiplier")
axes[1].legend()
plt.tight_layout()
plt.savefig("io_multipliers.png", dpi=150, bbox_inches="tight")
plt.close()
print("\nFigure saved: io_multipliers.png")
Step 3: Regional Convergence and Spatial Econometrics
import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
from scipy import stats
# -----------------------------------------------------------------# Simulate panel data: 50 regions, 20 years# -----------------------------------------------------------------
np.random.seed(123)
n_r = 50
T = 20# Initial per-capita income (log scale, heterogeneous)
ln_y0 = np.random.uniform(9.0, 11.5, n_r) # ln(initial income)# True beta: -0.05 (convergence), alpha: 0.5
beta_true = -0.05
alpha_true = 0.5
g = alpha_true + beta_true * ln_y0 + np.random.normal(0, 0.02, n_r) # annual growth rate# -----------------------------------------------------------------# Absolute beta-convergence OLS# -----------------------------------------------------------------
X_reg = sm.add_constant(ln_y0)
model = sm.OLS(g, X_reg)
result = model.fit(cov_type="HC3") # heteroskedasticity-robust SEprint("=== Absolute Beta-Convergence ===")
print(result.summary().tables[1])
beta_hat = result.params[1]
lam = -np.log(1 + beta_hat * T) / T if (1 + beta_hat * T) > 0else np.nan
half_life = np.log(2) / lam if lam > 0else np.inf
print(f"\nEstimated beta: {beta_hat:.4f}")
print(f"Convergence speed lambda: {lam:.4f}")
print(f"Half-life: {half_life:.1f} years")
# -----------------------------------------------------------------# Sigma-convergence: dispersion over time# -----------------------------------------------------------------
sigma = []
for t inrange(T):
ln_yt = ln_y0 + g * t # simplified linear approximation
sigma.append(ln_yt.std())
sigma = np.array(sigma)
# -----------------------------------------------------------------# Krugman Specialization Index between all region pairs# -----------------------------------------------------------------
n_regions_k = 8
n_industries_k = 6
emp_k = np.random.randint(500, 5000, (n_regions_k, n_industries_k)).astype(float)
shares = emp_k / emp_k.sum(axis=1, keepdims=True) # industry shares per region
K_matrix = np.zeros((n_regions_k, n_regions_k))
for i inrange(n_regions_k):
for j inrange(n_regions_k):
K_matrix[i, j] = 0.5 * np.abs(shares[i] - shares[j]).sum()
print(f"\nKrugman Specialization Index (mean): {K_matrix[np.triu_indices(n_regions_k, k=1)].mean():.3f}")
# -----------------------------------------------------------------# Economic base multiplier# -----------------------------------------------------------------# Using LQ > 1 rule to identify basic employment
LQ_k = location_quotient(emp_k) # reuse function from Step 1# Basic employment in each regiondefeconomic_base_multiplier(emp, lq):
"""Compute economic base multiplier for each region."""
total_emp = emp.sum(axis=1)
basic_emp = np.where(lq > 1.0, emp - emp.sum(axis=0) / emp.sum() * emp.sum(axis=1, keepdims=True), 0)
basic_emp = np.maximum(basic_emp, 0).sum(axis=1)
multiplier = total_emp / np.where(basic_emp > 0, basic_emp, 1)
return multiplier, basic_emp
mult_k, basic_k = economic_base_multiplier(emp_k, LQ_k)
for r, (m, b) inenumerate(zip(mult_k, basic_k)):
print(f"Region {chr(65+r)}: Basic emp = {b:.0f}, Multiplier = {m:.2f}")
# -----------------------------------------------------------------# Visualization# -----------------------------------------------------------------
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
# Beta-convergence scatter
axes[0].scatter(ln_y0, g, alpha=0.6, edgecolors="k", linewidths=0.4)
x_line = np.linspace(ln_y0.min(), ln_y0.max(), 100)
axes[0].plot(x_line, result.params[0] + result.params[1] * x_line, "r-", lw=2)
axes[0].set_xlabel("ln(Initial Income)")
axes[0].set_ylabel("Average Annual Growth Rate")
axes[0].set_title(f"β-Convergence (β={beta_hat:.3f})")
# Sigma-convergence
axes[1].plot(range(T), sigma, "o-", color="steelblue")
axes[1].set_xlabel("Year")
axes[1].set_ylabel("Std Dev of ln(Income)")
axes[1].set_title("σ-Convergence")
# Krugman specialization matrix
im2 = axes[2].imshow(K_matrix, cmap="YlOrRd", vmin=0, vmax=0.6)
axes[2].set_title("Krugman Specialization Index")
plt.colorbar(im2, ax=axes[2], label="K")
plt.tight_layout()
plt.savefig("regional_convergence.png", dpi=150, bbox_inches="tight")
plt.close()
print("\nFigure saved: regional_convergence.png")
deflocation_quotient(emp_matrix):
"""Reusable LQ function."""
total_regional = emp_matrix.sum(axis=1, keepdims=True)
total_national = emp_matrix.sum()
industry_national = emp_matrix.sum(axis=0, keepdims=True)
share_regional = emp_matrix / total_regional
share_national = industry_national / total_national
return share_regional / share_national
Advanced Usage
Esteban-Ray Polarization Index
import numpy as np
defesteban_ray_polarization(income, population, alpha=1.6):
"""Compute Esteban-Ray (1994) polarization index.
Args:
income: array of group mean incomes
population: array of group population shares
alpha: polarization sensitivity parameter (1 <= alpha <= 1.6)
Returns:
polarization index P
"""
n = len(income)
P = 0.0for i inrange(n):
for j inrange(n):
P += population[i]**(1 + alpha) * population[j] * abs(income[i] - income[j])
return P
# Example: 5 income groups
income_groups = np.array([15000, 30000, 50000, 80000, 150000])
pop_shares = np.array([0.20, 0.25, 0.30, 0.15, 0.10])
P = esteban_ray_polarization(income_groups, pop_shares)
print(f"Esteban-Ray Polarization: {P:.2f}")
Spatial Autocorrelation of Regional Income
import numpy as np
from scipy.spatial.distance import cdist
defmorans_i_regional(values, coords, k_neighbors=5):
"""Moran's I with k-nearest-neighbor spatial weights.
Args:
values: (n,) array of regional values
coords: (n, 2) array of region centroids
k_neighbors: number of neighbors
Returns:
morans_i, expected_i, z_score
"""
n = len(values)
dist_matrix = cdist(coords, coords)
np.fill_diagonal(dist_matrix, np.inf)
# Build binary KNN weight matrix
W = np.zeros((n, n))
for i inrange(n):
knn_idx = np.argsort(dist_matrix[i])[:k_neighbors]
W[i, knn_idx] = 1
W_sum = W.sum()
# Moran's I
z = values - values.mean()
numer = n * (W * np.outer(z, z)).sum()
denom = W_sum * (z**2).sum()
I = numer / denom
E_I = -1 / (n - 1)
# Variance under normality assumption
S1 = 0.5 * ((W + W.T)**2).sum()
S2 = ((W.sum(axis=1) + W.sum(axis=0))**2).sum()
n2 = n**2
num_var = n * ((n2 - 3*n + 3)*S1 - n*S2 + 3*W_sum**2)
den_var = (n-1)*(n-2)*(n-3)*W_sum**2
kurtosis = (z**4).mean() / (z**2).mean()**2
var_I = num_var / den_var - kurtosis * ((n2-n)*S1 - 2*n*S2 + 6*W_sum**2) / den_var
z_score = (I - E_I) / np.sqrt(abs(var_I))
return I, E_I, z_score
np.random.seed(42)
n_reg = 30
coords = np.random.uniform(0, 100, (n_reg, 2))
# Spatially autocorrelated income
income = 50000 + 10000 * np.sin(coords[:, 0] / 20) + np.random.normal(0, 2000, n_reg)
I_stat, E_I, z = morans_i_regional(income, coords)
print(f"Moran's I = {I_stat:.4f}, E[I] = {E_I:.4f}, Z = {z:.2f}")
Gravity Model for Interregional Trade
import numpy as np
import pandas as pd
import statsmodels.api as sm
defgravity_model(gdp_i, gdp_j, dist_ij, trade_ij):
"""Estimate gravity model: ln(T_ij) = a + b*ln(GDP_i) + c*ln(GDP_j) + d*ln(dist_ij).
Args:
gdp_i, gdp_j: (n,) arrays of exporter/importer GDP
dist_ij: (n,) distance between regions
trade_ij: (n,) bilateral trade flows
Returns:
OLS result
"""
df = pd.DataFrame({
"ln_trade": np.log(trade_ij),
"ln_gdp_i": np.log(gdp_i),
"ln_gdp_j": np.log(gdp_j),
"ln_dist": np.log(dist_ij),
}).dropna()
X = sm.add_constant(df[["ln_gdp_i", "ln_gdp_j", "ln_dist"]])
model = sm.OLS(df["ln_trade"], X)
return model.fit(cov_type="HC3")
# Simulate 200 region pairs
np.random.seed(42)
n_pairs = 200
gdp_i = np.random.lognormal(10, 1.5, n_pairs)
gdp_j = np.random.lognormal(10, 1.5, n_pairs)
dist_ij = np.random.uniform(50, 2000, n_pairs)
# True gravity: T = GDP_i^0.8 * GDP_j^0.9 * dist^-1.2 * noise
trade_ij = (gdp_i**0.8 * gdp_j**0.9 * dist_ij**(-1.2) *
np.exp(np.random.normal(0, 0.3, n_pairs)))
result = gravity_model(gdp_i, gdp_j, dist_ij, trade_ij)
print("Gravity Model Estimates:")
print(result.summary().tables[1])
Troubleshooting
Problem
Cause
Fix
np.linalg.inv singular matrix
Near-zero or zero column in A
Check for all-zero industries; use np.linalg.lstsq or pseudoinverse
Negative basic employment
LQ < 1 or rounding
np.maximum(basic_emp, 0) clamp
Convergence beta insignificant
Too few regions or time span
Increase sample; use club-convergence tests
Division by zero in LQ
Region has zero employment
Filter out empty regions before computing
IO multiplier > 10
Unrealistic A matrix
Check that column sums of A < 1 (no super-multiplier)
Moran's I close to -1/(n-1)
No spatial pattern
Expected; test for significance with permutation
External Resources
Isard, W. (1960). Methods of Regional Analysis. MIT Press.
Dunn, E. S. (1960). A statistical and analytical technique for regional analysis. Papers in Regional Science.
Leontief, W. (1986). Input-Output Economics. Oxford University Press.
Anselin, L. (1988). Spatial Econometrics. Kluwer Academic.
Esteban, J., & Ray, D. (1994). On the measurement of polarization. Econometrica, 62(4), 819-851.
