Skip to content

Model card — Functional shocks (FPCA, FLP, FVAR)

functional_pca · flp · flp_scenario · fvar_scenario

Most impulse-response analysis assumes the shock is a number — a 25bp surprise, a one-standard-deviation innovation. But some of the most important shocks are curves: a monetary announcement does not move "the" interest rate, it moves the whole yield curve, and two announcements with the same short-rate surprise can twist the curve in entirely different ways. The Inoue-Rossi (2021) idea implemented here: treat the shock as a function, compress the observed panel of curves into a few principal-component scores, estimate how the outcome responds to each score, and then ask about any scenario — "the whole curve shifts up 100bp", "the curve flattens" — by projecting that scenario onto the score space.

The pipeline is

curves (T x M)  --FPCA-->  scores (T x K)  --LP or VAR-->  response of y
                                                            to a curve scenario

and its trust anchor is a reconstruction identity you can verify in one line: feeding the j-th eigenfunction back in as the scenario returns exactly the j-th score's coefficient path.

functional_pca — principal components of a curve panel

What it estimates. The functional principal-component decomposition of a T × M panel of curves (each row one period's curve on a common M-point grid): the cross-sectional mean_curve, the top n_factors eigenfunctions (orthonormal K × M loadings), the scores (T × K — each period's curve expressed in eigenfunction coordinates, demeaned), eigenvalues, and the explained variance shares. For yield curves the first two components are reliably "level" and "slope" — the empirical basis the rest of the family works in.

Assumptions. Curves observed on the same grid every period, no missing points; enough factors that the retained components capture the variation you care about (check explained). PCA is a linear decomposition — a curve panel whose variation is genuinely nonlinear in maturity will smear across many components.

When to use (and when not). Use it as the first stage of the functional- shock pipeline, or standalone as a disciplined dimension reduction for any curve-valued series (yield curves, forward curves, distributions binned to a grid). Not for panels of unrelated series (use factor models on standardized data), and not when M changes over time — interpolate to a fixed grid first.

Key arguments and defaults. n_factors = 3 — for yield curves the level/slope/curvature convention; raise it until explained sums close enough to 1 for your purposes. Sign convention (documented and stable): each eigenfunction is normalized so its largest-magnitude entry is positive — eigenvectors are only defined up to sign, and pinning it makes scores and downstream betas reproducible across runs and platforms.

How to read the output. eigenfunctions[j] is the shape of the j-th mode of variation across the grid (all-positive ≈ level; monotone sign-switching ≈ slope); scores[t, j] is how much of that shape period t's curve contains (relative to the mean curve); explained[j] is the fraction of total variance it accounts for. Reconstruction: curves[t] ≈ mean_curve + scores[t] @ eigenfunctions.

Failure modes. Too few factors silently truncates the shock space — a scenario orthogonal to the retained eigenfunctions has weights near zero and a spuriously flat "response" (see flp_scenario). Wildly different units across grid points let one region of the curve dominate the decomposition. Near-equal eigenvalues make individual eigenfunctions unstable (their span is fine, their labels are not).

Validated against. numpy.linalg.eigh on the sample covariance of the demeaned panel at 1e-8 — eigenvalues, explained shares, and (after the documented sign normalization) eigenfunctions and scores (fixtures/tsecon-funcshock.json, generated by fixtures/generate_tsecon-funcshock_fixtures.py).

References. Inoue & Rossi (2021, Quantitative Economics 12:1085-1138); Ramsay & Silverman (2005, Functional Data Analysis).

flp — functional local projections

What it estimates. Per horizon h, a joint regression of y_{t+h} on all K scores at once, plus a constant and n_lag_controls lags of y:

y_{t+h} = beta_h' scores_t + c_h + sum_l phi_{h,l} y_{t-l} + u_{t,h}

with Newey-West HAC standard errors. betas[h] is the length-K vector of score coefficients at horizon h; jointness matters because the full K × K covariance covs[h] is what lets flp_scenario put a standard error on any linear combination of the score paths.

Assumptions. The scores are valid impulse variables — exogenous conditional on the lag controls. That is a real identification assumption, not a technicality: raw yield-curve scores are equilibrium objects. The clean Inoue-Rossi use case feeds curves of changes around identified events (monetary announcement windows) so the scores inherit the event-study identification. Serial correlation in the overlapping residuals is handled by HAC, as in any LP.

When to use (and when not). Use it when you already have scores (from functional_pca or elsewhere) and want the per-score response paths and their joint covariance. If your question is "how does y respond to this curve scenario?", call flp_scenario directly — it runs FPCA and FLP internally and does the projection arithmetic for you.

Key arguments and defaults. horizons = 8, n_lag_controls = 2; hac_maxlags defaults to h + n_lag_controls per horizon (bandwidth grows with the overlap, the standard LP-HAC prescription).

How to read the output. betas is (H+1) × K (row h, column k = the h-step response of y to a one-unit innovation in score k); se its per-element standard errors; covs[h] the joint K × K covariance at horizon h; nobs the per-horizon sample sizes. Score units are curve units times the eigenfunction normalization — economically meaningful statements come from flp_scenario's weighted combinations, not from a raw score coefficient.

Failure modes. Endogenous scores (the identification point above). Near-collinear scores — by FPCA construction in-sample scores are orthogonal, but scores from sub-samples or externally supplied ones need not be. Too many factors relative to T inflates every horizon's regression.

Validated against. statsmodels OLS with Newey-West HAC on the identical joint design, horizon by horizon (fixtures/tsecon-funcshock.json).

References. Inoue & Rossi (2021); Jordà (2005); Newey & West (1987).

flp_scenario — response of y to a whole-curve scenario

What it estimates. The impulse response of y to the entire curve moving by delta (a length-M vector — your scenario, in curve units). It runs functional_pca(curves, n_factors), projects the scenario onto the eigenfunctions to get weights w = Phi @ delta, runs the joint flp, and returns

response_h = w' beta_h,        se_h = sqrt( w' covs[h] w )

— the response to the scenario as a linear combination of score responses, with the correct joint-covariance standard error.

The reconstruction identity — the trust anchor. Because the eigenfunctions are orthonormal, setting delta equal to the j-th eigenfunction makes the weights exactly the j-th unit vector, so the scenario response must equal the j-th column of flp's betas to machine precision. This is pinned in the test suite and verifiable in one line (shown in the example below): if the identity holds, the FPCA projection, the joint regression, and the scenario arithmetic are all wired together correctly.

When to use (and when not). This is the headline function of the family: "what happens to output if the whole yield curve shifts up 100bp?", "…if it flattens the way it did around this announcement?". Use delta vectors that are economically named (parallel shift, observed announcement-window change, an eigenfunction). Do not use scenarios far outside the span of the retained eigenfunctions — check weights against explained: a scenario the factor space cannot represent gets silently truncated to its projection.

Key arguments and defaults. As functional_pca + flp: n_factors = 3, horizons = 8, n_lag_controls = 2, hac_maxlags = h + p default. delta is in the units of the curves (e.g. percentage points of yield at each maturity).

How to read the output. response[h] ± z·se[h] is the IRF of y to the scenario; weights tells you how the scenario decomposes across eigenfunctions (and how much of it was representable); betas and explained are passed through for inspection.

Failure modes. A scenario orthogonal to the factor space (tiny weights despite a large delta — raise n_factors or admit the data never exhibits that shape). Everything from flp (score endogeneity above all). And scale errors: delta in basis points when curves are in percent produces responses 100× too large.

Validated against. The reconstruction identity (scenario = j-th eigenfunction ⇒ response = j-th beta path, machine precision) is test-pinned, on top of the functional_pca and flp validations above.

References. Inoue & Rossi (2021, Quantitative Economics 12:1085-1138).

fvar_scenario — the VAR route to the same question

What it estimates. The same whole-curve scenario response, but through a VAR instead of local projections: fit a VAR(lags) on [scores_1..K, y] with the scores ordered first, identify by Cholesky, set the score-block structural innovation to w = Phi' delta with the outcome's own structural shock zero, and iterate the VAR forward.

The identification caveat, plainly. Cholesky with scores first means: the scenario is allowed to move y on impact, but y's own shock moves no score on impact. The impact response of y (implied_outcome_innovation — the Cholesky-implied comovement) is therefore a modeling assumption, not an estimate robust to ordering. If the outcome in truth affects the curve within the period (an activity indicator moving yields intraday), the FVAR impact response is contaminated. flp_scenario makes the analogous assumption through its regression timing but is robust to misspecification of the long-run dynamics; agreement between the two routes is the useful diagnostic, not either one alone.

When to use (and when not). Use it as the cross-check on flp_scenario (close agreement is comforting; divergence at long horizons usually indicts the VAR's extrapolation, at short horizons the ordering), or when you want the VAR's tight parametric long-horizon behavior from a short sample. Do not present it alone as identification-free — see the caveat.

Key arguments and defaults. n_factors = 3, lags = 2, horizon = 10. No standard errors are returned on this route — it is a point-estimate cross-check.

How to read the output. response_outcome[h] is the scenario IRF of y; responses is (H+1) × (K+1) with the score paths first and the outcome last, so you can also see how the curve itself propagates; implied_outcome_innovation is the impact comovement discussed above (= response_outcome[0]); weights as in flp_scenario.

Failure modes. The ordering assumption (above). Unstable VAR estimates on short samples propagate errors into every horizon — unlike LP, a bad one-step-ahead model poisons the whole path. Too few lags leave residual autocorrelation that Cholesky mistakes for contemporaneous structure.

Validated against. statsmodels VAR([scores, y]).fit(lags, trend="c") plus orth_ma_rep (Cholesky-orthogonalized MA representation) and a scipy triangular solve, at 1e-8 (fixtures/tsecon-funcshock.json); the impact identity responses[0][:K] == weights is property-tested.

References. Inoue & Rossi (2021); Sims (1980).

Worked example — a yield-curve panel end to end

Simulated Nelson-Siegel yields at 8 maturities driven by AR(1) level and slope factors; activity y falls when the level of the curve rises (y_t = 0.6 y_{t-1} - 0.4 level_t + noise):

import numpy as np, tsecon

rng = np.random.default_rng(0)
n = 400
mats = np.array([0.25, 0.5, 1, 2, 3, 5, 7, 10])       # maturities (years), M = 8
lam = 0.7
slope_load = (1 - np.exp(-lam * mats)) / (lam * mats)  # Nelson-Siegel slope loading

level = np.zeros(n); slope = np.zeros(n)
ul = rng.standard_normal(n); us = rng.standard_normal(n)
for t in range(1, n):
    level[t] = 0.9 * level[t - 1] + 0.3 * ul[t]
    slope[t] = 0.7 * slope[t - 1] + 0.3 * us[t]
curves = (4.0 + np.outer(level, np.ones(len(mats)))
          + np.outer(slope, slope_load)
          + 0.05 * rng.standard_normal((n, len(mats))))   # T x M yield panel

# activity falls when the LEVEL of the curve rises
y = np.zeros(n); ey = rng.standard_normal(n)
for t in range(1, n):
    y[t] = 0.6 * y[t - 1] - 0.4 * level[t] + 0.2 * ey[t]

# ---- 1. functional_pca: level and slope recovered from the panel
p = tsecon.functional_pca(curves, n_factors=2)
print("explained:", np.round(np.asarray(p["explained"]), 4))
phi = np.asarray(p["eigenfunctions"])                  # K x M
print("phi_1 (level-like):", np.round(phi[0], 3))
print("phi_2 (slope-like):", np.round(phi[1], 3))
scores = np.asarray(p["scores"])                       # T x K
print("corr(score_1, true level):", round(float(np.corrcoef(scores[:, 0], level)[0, 1]), 3))

# ---- 2. flp: joint per-horizon response of y to the scores
r = tsecon.flp(y, scores, horizons=8, n_lag_controls=2)
betas = np.asarray(r["betas"]); se = np.asarray(r["se"])   # (H+1) x K
print("\nbeta on score_1, h=0..4:", np.round(betas[:5, 0], 3))
print("se,               h=0..4:", np.round(se[:5, 0], 3))

# ---- 3. the reconstruction identity (trust anchor)
s1 = tsecon.flp_scenario(y, curves, phi[0], n_factors=2, horizons=8, n_lag_controls=2)
print("\nscenario = phi_1  =>  weights:", np.round(np.asarray(s1["weights"]), 6))
print("max |response - beta path 1| :",
      float(np.max(np.abs(np.asarray(s1["response"]) - betas[:, 0]))))

# ---- 4. an economically named scenario: a 100bp parallel shift
shift = np.ones(len(mats))                              # +1.0 at every maturity
s = tsecon.flp_scenario(y, curves, shift, n_factors=2, horizons=8, n_lag_controls=2)
print("\nparallel-shift weights:", np.round(np.asarray(s["weights"]), 3))
print("response h=0..8:", np.round(np.asarray(s["response"]), 3))
print("se       h=0..8:", np.round(np.asarray(s["se"]), 3))

# ---- 5. the FVAR route on the same scenario
v = tsecon.fvar_scenario(y, curves, shift, n_factors=2, lags=2, horizon=8)
print("\nFVAR response :", np.round(np.asarray(v["response_outcome"]), 3))
print("implied outcome innovation (Cholesky comovement):",
      round(float(v["implied_outcome_innovation"]), 3))
# explained: [0.9738 0.0224]
# phi_1 (level-like): [0.389 0.382 0.373 0.356 0.345 0.329 0.326 0.321]
# phi_2 (slope-like): [ 0.499  0.394  0.253  0.028 -0.152 -0.325 -0.424 -0.472]
# corr(score_1, true level): 0.949
#
# beta on score_1, h=0..4: [-0.134 -0.191 -0.214 -0.206 -0.183]
# se,               h=0..4: [0.007 0.013 0.016 0.02  0.024]
#
# scenario = phi_1  =>  weights: [1. 0.]
# max |response - beta path 1| : 5.551115123125783e-17
#
# parallel-shift weights: [ 2.821 -0.199]
# response h=0..8: [-0.413 -0.597 -0.673 -0.642 -0.56  -0.512 -0.481 -0.426 -0.361]
# se       h=0..8: [0.024 0.038 0.047 0.058 0.071 0.077 0.076 0.074 0.083]
#
# FVAR response : [-0.416 -0.597 -0.659 -0.662 -0.631 -0.583 -0.527 -0.469 -0.413]
# implied outcome innovation (Cholesky comovement): -0.416

Everything a user should check is on display: FPCA recovers a flat, all-positive level eigenfunction and a sign-switching slope one (97.4% and 2.2% of variance); the reconstruction identity holds to 6e-17; the 100bp parallel shift — weights (2.82, −0.20), almost pure level — depresses activity with the true hump shape (impact ≈ −0.41, matching the DGP's −0.4, peak ≈ −0.67 at h = 2); and the FVAR route, with its scores-first Cholesky assumption, tracks the FLP route closely at every horizon — the agreement diagnostic doing its job.