/20-wenddymacro-python-econ-skill

name: python-econ-computing
description: Use when writing Python code for DSGE models, HANK models, numerical economic computation, causal inference, or quantitative economic data analysis

Python Economic Numerical Computing

• Author：Wenli Xu
• Email： wlxu@cityu.edu.mo
• 2026-03-11

---

Overview

Best practices for macroeconomic modeling (DSGE/HANK), causal inference, and data analysis in Python. Core principle: **vectorize first, accelerate loops with Numba, keep code structure aligned with economic theory**.

---

Library Quick Reference

| Use Case | Preferred Libraries | |----------|-------------------| | Numerical core | `numpy`, `scipy` | | Loop acceleration | `numba` (`@njit`, `@njit(parallel=True)`) | | Economics toolkit | `quantecon` | | HANK / sequence space | `sequence_jacobian` (SSJ) | | Heterogeneous agents | `HARK` | | **Linear models with FE** | **`pyfixest`** (`pip install pyfixest`) | | **DID / DD / DDD** | **`diff-diff`** (`pip install diff-diff`) | | **IV / 2SLS / GMM** | **`linearmodels`** (or `pyfixest` for panel IV with FE) | | **RD / RDD / RKD** | **`rdrobust`**, `rddensity`, `rdlocrand` | | **Synthetic Control** | **`pysynth`**, `synth_control`, `sdid` | | **Matching** | **`causalml`**, `pymatch`, `econml` | | **Causal ML / DML** | **`econml`**, `dowhy` | | Data manipulation | `pandas`, `polars` (large datasets) | | Visualization | `matplotlib`, `seaborn` |

---

DSGE Models

Linearization and Solution (Blanchard-Kahn)

import numpy as np
from scipy.linalg import ordqz

def solve_bk(A, B, n_fwd):
    """
    Solve linear DSGE: A E_t[x_{t+1}] = B x_t + C eps_t
    n_fwd: number of forward-looking variables
    Returns decision rule matrix P such that x_t = P x_{t-1} + ...
    """
    AA, BB, alpha, beta, Q, Z = ordqz(A, B, sort='ouc')
    n = A.shape[0]
    Z21 = Z[n - n_fwd:, :n - n_fwd]
    Z22 = Z[n - n_fwd:, n - n_fwd:]
    P = -np.linalg.solve(Z22, Z21)
    return P

Perturbation Methods (Second-Order Approximation)

• Use `quantecon.lqcontrol` for LQ problems
• Higher-order perturbation: `perturbpy` or manual implementation
• Steady-state solving: `scipy.optimize.fsolve` / `root`

---

HANK Models

Sequence-Space Jacobian Method (SSJ)

import sequence_jacobian as sj

# 1. Define steady-state blocks
@sj.simple
def household_ss(r, w, beta, sigma):
    # Return steady-state aggregates
    ...

# 2. Build DAG
model = sj.create_model([household_block, firm_block, market_clearing],
                         name='HANK')

# 3. Solve steady state
ss = model.solve_steady_state(calibration, unknowns, targets)

# 4. Compute Jacobians → solve transition dynamics
G = model.solve_jacobian(ss, unknowns, targets, T=300)

Value Function Iteration — Numba Accelerated

from numba import njit
import numpy as np

@njit
def vfi(V0, a_grid, y_grid, r, beta, sigma, tol=1e-8, max_iter=1000):
    """Heterogeneous agent VFI over asset grid × income grid"""
    n_a, n_y = len(a_grid), len(y_grid)
    V = V0.copy()
    policy = np.zeros((n_a, n_y))

    for it in range(max_iter):
        V_new = np.empty_like(V)
        for ia in range(n_a):
            for iy in range(n_y):
                best_val = -1e10
                best_a = 0
                for ia2 in range(n_a):
                    c = (1 + r) * a_grid[ia] + y_grid[iy] - a_grid[ia2]
                    if c <= 0:
                        continue
                    u = c ** (1 - sigma) / (1 - sigma)
                    val = u + beta * V[:, iy].mean()  # use transition matrix in practice
                    if val > best_val:
                        best_val = val
                        best_a = ia2
                V_new[ia, iy] = best_val
                policy[ia, iy] = a_grid[best_a]
        if np.max(np.abs(V_new - V)) < tol:
            break
        V = V_new
    return V, policy

Distribution Iteration (Young 2010)

def iterate_distribution(policy_idx, trans_mat, dist0, T=500):
    """Iterate joint distribution to steady state given policy indices and income transition matrix"""
    dist = dist0.copy()
    n_a, n_y = dist.shape
    for _ in range(T):
        dist_new = np.zeros_like(dist)
        for iy in range(n_y):
            for iy2 in range(n_y):
                dist_new[policy_idx[:, iy], iy2] += dist[:, iy] * trans_mat[iy, iy2]
        dist = dist_new
    return dist

---

Linear Models with Fixed Effects (pyfixest)

**Rule: For any OLS/Poisson/Logit with fixed effects, use `pyfixest`. It mirrors R's `fixest` syntax.**

import pyfixest as pf

# OLS with unit + time FE, cluster-robust SEs
fit = pf.feols("y ~ treat_post | unit + year",
               data=df, vcov={"CRV1": "id"})
fit.summary()

# Multiple high-dimensional FE (Frisch-Waugh absorbed)
fit = pf.feols("y ~ x1 + x2 | unit + year + industry",
               data=df, vcov={"CRV1": "id"})

# Wild cluster bootstrap (few clusters, <50)
fit = pf.feols("y ~ treat_post | unit + year",
               data=df, vcov={"CRV1": "id"})
fit.wildboottest(param="treat_post", B=9999, seed=42)

# Event study via i() syntax
fit = pf.feols("y ~ i(rel_year, ref=-1) | unit + year",
               data=df, vcov={"CRV1": "id"})
pf.iplot(fit)  # event study plot

# Poisson (count / log-linear) with FE
fit_pois = pf.fepois("y ~ treat_post | unit + year",
                     data=df, vcov={"CRV1": "id"})

# Access results
fit.coef()           # coefficient estimates
fit.se()             # standard errors
fit.pvalue()         # p-values
fit.confint()        # confidence intervals
fit._N               # number of observations

pyfixest vs statsmodels

| Use case | Use | |----------|-----| | OLS / WLS with any FE | `pyfixest` | | Poisson / logit with FE | `pyfixest` | | Wild bootstrap | `pyfixest` | | Time-series ARIMA, VAR | `statsmodels` | | MLE / GLM without FE | `statsmodels` |

---

Causal Inference: DID / DD / DDD Methods

**Rule: For any DiD, DD, DDD, or s