xjtulyc Regional Economics Analysis
When to Use This Skill
Use this skill when you need to:
- Analyze regional economic structure and competitive advantage
- Decompose employment/output changes into national, industry, and competitive effects
- Measure industrial specialization and diversification
- Build regional input-output tables and compute multipliers
- Estimate spatial econometric models for regional data
- Assess regional convergence and divergence patterns
Trigger keywords: shift-share analysis, location quotient, economic base theory, input-output model, regional multiplier, Krugman specialization, beta convergence, sigma convergence, spatial lag model, spatial error model, regional growth, industrial cluster, agglomeration economies, regional GDP, interregional trade.
Background & Key Concepts
Location Quotient
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:
$$\Delta e_{ir} = \underbrace{e_{ir} \cdot g_n}{\text{National Share}} + \underbrace{e{ir}(g_{in} - g_n)}{\text{Industry Mix}} + \underbrace{e{ir}(g_{ir} - g_{in})}_{\text{Competitive Effect}}$$
where $g_n$, $g_{in}$, $g_{ir}$ are national, national-industry, and regional-industry growth rates.
Economic Base Multiplier
Basic sector employment $B$ drives total employment $T$ via multiplier $k$:
$$k = \frac{T}{B}, \quad B = \sum_i \max!\left(e_{ir} - \frac{e_{in}}{E_n} E_r,\ 0\right)$$
Input-Output Leontief Model
Given technical coefficient matrix $\mathbf{A}$ (where $a_{ij} = z_{ij}/X_j$):
$$\mathbf{X} = (\mathbf{I} - \mathbf{A})^{-1} \mathbf{f}$$
The Leontief inverse $\mathbf{L} = (\mathbf{I}-\mathbf{A})^{-1}$ gives total output multipliers. Output multiplier for sector $j$: $m_j = \sum_i L_{ij}$.
Krugman Specialization Index
$$K_{rs} = \frac{1}{2}\sum_i \left|\frac{e_{ir}}{E_r} - \frac{e_{is}}{E_s}\right| \in [0,1]$$
$K=0$: identical structure; $K=1$: completely different.
Beta-Convergence
Absolute $\beta$-convergence regression:
$$\frac{1}{T}\ln\frac{y_{r,T}}{y_{r,0}} = \alpha + \beta \ln y_{r,0} + \varepsilon_r$$
$\beta < 0$ indicates poorer regions grow faster (convergence). Implied convergence speed: $\lambda = -\ln(1+\beta)/T$.
Environment Setup
pip install numpy>=1.24 pandas>=2.0 scipy>=1.11 statsmodels>=0.14 \
matplotlib>=3.7 geopandas>=0.14
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
# -----------------------------------------------------------------
# Simulate regional employment data
# 8 regions, 6 industries, 2 time periods (t0, t1)
# -----------------------------------------------------------------
np.random.seed(42)
n_regions = 8
n_industries = 6
region_names = [f"Region_{r}" for r in "ABCDEFGH"]
industry_names = ["Manufacturing", "Finance", "Retail",
"Construction", "Healthcare", "Tech"]
# Base employment (t0)
emp_t0 = np.random.randint(500, 5000, (n_regions, n_industries)).astype(float)
# National industry growth rates (differential by industry)
nat_industry_growth = np.array([0.02, 0.08, 0.03, 0.04, 0.06, 0.12])
# Regional competitive advantages (slight variation around national)
comp_effect = np.random.uniform(-0.03, 0.05, (n_regions, n_industries))
# t1 employment
emp_t1 = emp_t0 * (1 + nat_industry_growth + comp_effect)
# -----------------------------------------------------------------
# Location Quotient
# -----------------------------------------------------------------
def location_quotient(emp_matrix):
"""Compute LQ for each region-industry cell.
Args:
emp_matrix: (n_regions, n_industries) array
Returns:
lq: same-shape array of location quotients
"""
total_regional = emp_matrix.sum(axis=1, keepdims=True) # E_r
total_national = emp_matrix.sum() # E_n
industry_national = emp_matrix.sum(axis=0, keepdims=True) # e_in
share_regional = emp_matrix / total_regional # e_ir/E_r
share_national = industry_national / total_national # e_in/E_n
lq = share_regional / share_national
return lq
lq = location_quotient(emp_t0)
lq_df = pd.DataFrame(lq, index=region_names, columns=industry_names)
print("Location Quotients (t0):")
print(lq_df.round(2))
print(f"\nSpecialized cells (LQ>1.25): {(lq > 1.25).sum()}")
# -----------------------------------------------------------------
# Classic Shift-Share (Dunn 1960)
# -----------------------------------------------------------------
def shift_share(emp_t0, emp_t1):
"""Dunn shift-share decomposition.
Returns national_share, industry_mix, competitive_effect
all shape (n_regions, n_industries).
"""
E_n0 = emp_t0.sum()
E_n1 = emp_t1.sum()
g_n = (E_n1 - E_n0) / E_n0 # national growth rate
e_in0 = emp_t0.sum(axis=0) # national industry emp t0
e_in1 = emp_t1.sum(axis=0) # national industry emp t1
g_in = (e_in1 - e_in0) / e_in0 # national industry growth rates
g_ir = (emp_t1 - emp_t0) / emp_t0 # regional-industry growth
ns = emp_t0 * g_n # national share effect
im = emp_t0 * (g_in - g_n) # industry mix effect
ce = emp_t0 * (g_ir - g_in) # competitive effect
return ns, im, ce
ns, im, ce = shift_share(emp_t0, emp_t1)
# Aggregate by region
ss_df = pd.DataFrame({
"National_Share": ns.sum(axis=1),
"Industry_Mix": im.sum(axis=1),
"Competitive": ce.sum(axis=1),
"Total_Change": (emp_t1 - emp_t0).sum(axis=1)
}, index=region_names)
print("\nShift-Share Decomposition (aggregate by region):")
print(ss_df.round(0))
# -----------------------------------------------------------------
# Visualization: Stacked bar chart
# -----------------------------------------------------------------
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# LQ heatmap
im_ax = axes[0].imshow(lq, cmap="RdYlGn", vmin=0.5, vmax=2.0, aspect="auto")
axes[0].set_xticks(range(n_industries)); axes[0].set_xticklabels(industry_names, rotation=30, ha="right")
axes[0].set_yticks(range(n_regions)); axes[0].set_yticklabels(region_names)
axes[0].set_title("Location Quotients")
plt.colorbar(im_ax, ax=axes[0], label="LQ")
# Shift-share stacked bar
x = np.arange(n_regions)
w = 0.6
axes[1].bar(x, ss_df["National_Share"], w, label="National Share", color="steelblue")
axes[1].bar(x, ss_df["Industry_Mix"], w, bottom=ss_df["National_Share"],
label="Industry Mix", color="orange")
axes[1].bar(x, ss_df["Competitive"], w,
bottom=ss_df["National_Share"] + ss_df["Industry_Mix"],
label="Competitive", color="green")
axes[1].axhline(0, color="black", linewidth=0.8)
axes[1].set_xticks(x); axes[1].set_xticklabels(region_names, rotation=30, ha="right")
axes[1].set_title("Shift-Share Decomposition")
axes[1].set_ylabel("Employment Change")
axes[1].legend()
plt.tight_layout()
plt.savefig("shift_share_analysis.png", dpi=150, bbox_inches="tight")
plt.close()
print("\nFigure saved: shift_share_analysis.png")
Step 2: Input-Output Model and Multipliers
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_j
print("\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 in zip(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 requirements
print("\nEmployment Multipliers:")
for s, m in zip(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 in zip(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 cells
for i in range(n):
for j in range(n):
axes[0].text(j, i, f"{L[i,j]:.2f}", ha="center", va="center",
color="white" if L[i,j] > 1.5 else "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 SE
print("=== 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) > 0 else np.nan
half_life = np.log(2) / lam if lam > 0 else 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 in range(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 in range(n_regions_k):
for j in range(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 region
def economic_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) in enumerate(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")
def location_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
def esteban_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.0
for i in range(n):
for j in range(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
def morans_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 in range(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
def gravity_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.
- PySAL documentation — Python library for spatial analysis
- FRED Regional Data — Federal Reserve regional economic data
- BEA Input-Output Tables
Examples
Example 1: Full Shift-Share Report for a State
import numpy as np
import pandas as pd
np.random.seed(2024)
n_industries = 10
industries = [
"Agriculture", "Mining", "Utilities", "Construction",
"Manufacturing", "Wholesale", "Retail", "Finance",
"Professional Services", "Healthcare"
]
# State employment vs. national employment
state_t0 = np.array([1200, 800, 300, 2500, 8000, 1800, 5000, 3200, 4500, 6000], dtype=float)
state_t1 = state_t0 * np.array([0.98, 0.85, 1.02, 1.05, 0.92, 1.08, 1.03, 1.12, 1.18, 1.15])
nat_t0 = state_t0 * np.random.uniform(15, 25, n_industries)
nat_t1 = nat_t0 * np.array([0.96, 0.80, 1.03, 1.06, 0.90, 1.10, 1.04, 1.15, 1.20, 1.13])
# Compute components
E_n0, E_n1 = nat_t0.sum(), nat_t1.sum()
g_n = (E_n1 - E_n0) / E_n0
g_in = (nat_t1 - nat_t0) / nat_t0
g_ir = (state_t1 - state_t0) / state_t0
report = pd.DataFrame({
"Industry": industries,
"Emp_t0": state_t0,
"Emp_t1": state_t1,
"LQ": (state_t0 / state_t0.sum()) / (nat_t0 / nat_t0.sum()),
"National_Share": state_t0 * g_n,
"Industry_Mix": state_t0 * (g_in - g_n),
"Competitive": state_t0 * (g_ir - g_in),
"Total_Change": state_t1 - state_t0,
})
report["Growing_Basic"] = (report["LQ"] > 1.25) & (report["Competitive"] > 0)
pd.set_option("display.float_format", "{:.0f}".format)
print("=== State Shift-Share Report ===")
print(report.to_string(index=False))
print(f"\nTotal employment change: {(state_t1 - state_t0).sum():.0f}")
print(f"Competitive effect total: {report['Competitive'].sum():.0f}")
Example 2: Input-Output Impact Assessment for Infrastructure Project
import numpy as np
import pandas as pd
# 6-sector IO table for a metropolitan region
sectors = ["Agriculture", "Manufacturing", "Construction",
"Trade", "Finance", "Services"]
n = len(sectors)
# Technical coefficient matrix (calibrated to reasonable values)
A = np.array([
[0.05, 0.08, 0.02, 0.01, 0.00, 0.01],
[0.10, 0.20, 0.15, 0.05, 0.02, 0.03],
[0.02, 0.05, 0.10, 0.03, 0.01, 0.02],
[0.08, 0.10, 0.08, 0.15, 0.05, 0.10],
[0.03, 0.05, 0.04, 0.06, 0.10, 0.08],
[0.12, 0.15, 0.10, 0.18, 0.20, 0.20],
])
# Verify sum < 1 (valid IO table)
assert all(A.sum(axis=0) < 1), "Column sums must be < 1"
L = np.linalg.inv(np.eye(n) - A)
# Infrastructure project: $500M construction spending
# Direct final demand shock
delta_f = np.zeros(n)
delta_f[2] = 500 # Construction sector
# Total output impact
delta_X = L @ delta_f
print("=== Infrastructure Project ($500M Construction) ===")
impact_df = pd.DataFrame({
"Sector": sectors,
"Direct_Demand": delta_f,
"Total_Output_Impact": delta_X,
"Multiplier": L[:, 2], # Construction column
})
print(impact_df.to_string(index=False))
print(f"\nTotal output multiplier for Construction: {L[:, 2].sum():.3f}")
print(f"Total economic impact: ${delta_X.sum():.1f}M")
print(f"Indirect/induced effects: ${delta_X.sum() - 500:.1f}M ({(delta_X.sum()/500-1)*100:.1f}%)")