diff --git a/docs/explanation/conditional-forecasting.md b/docs/explanation/conditional-forecasting.md new file mode 100644 index 0000000..a4523ae --- /dev/null +++ b/docs/explanation/conditional-forecasting.md @@ -0,0 +1,53 @@ +# Conditional Forecasting + +## The problem + +A standard VAR forecast projects all variables forward simultaneously, with no external constraints. But many questions in macroeconomics and finance require **conditional** projections: + +- Central banks publish forecasts conditioned on assumed interest rate paths +- Financial institutions stress-test portfolios under assumed macroeconomic scenarios +- Policy analysts ask "what if" questions about specific variable trajectories + +Conditional forecasting provides the mathematical framework for these exercises. + +## The idea + +A VAR forecast can be decomposed into two parts: + +$$y_{t+h} = \underbrace{y_{t+h}^{u}}_{\text{unconditional}} + \underbrace{\sum_{s=0}^{h} \Phi_{h-s} \, P \, \varepsilon_{s}}_{\text{shock-driven deviation}}$$ + +where $y^{u}_{t+h}$ is the unconditional forecast (no future shocks), $\Phi_j$ are the moving-average (MA) coefficient matrices, $P$ is the impact matrix (Cholesky factor of $\Sigma$ or a structural matrix), and $\varepsilon_s$ are the future structural shocks. + +The unconditional forecast is deterministic given the posterior draw. The future shocks are unknown. Conditional forecasting amounts to **choosing the shock paths** that make the forecast satisfy the desired constraints. + +## The Waggoner-Zha algorithm + +Waggoner & Zha (1999) showed that this can be formulated as a linear system. Stack all future shocks into a vector $\varepsilon = [\varepsilon_0, \varepsilon_1, \ldots, \varepsilon_{H-1}]$ of length $H \times n$, where $H$ is the number of forecast steps and $n$ is the number of variables. + +Each constraint (e.g., "variable $i$ equals value $v$ at period $h$") translates into a linear equation: + +$$R \, \varepsilon = c$$ + +where $R$ is constructed from the MA coefficients and the impact matrix, and $c$ contains the differences between target values and unconditional forecasts. + +The minimum-norm solution $\varepsilon^* = R^\top (R R^\top)^{-1} c$ gives the **smallest set of shocks** (in a least-squares sense) that satisfies all constraints. This is computed via `numpy.linalg.lstsq`. + +The conditional forecast is then: + +$$y_{t+h}^{c} = y_{t+h}^{u} + \sum_{s=0}^{h} \Phi_{h-s} \, P \, \varepsilon^*_s$$ + +## Structural extensions + +When the model is identified (you have a structural impact matrix $P$ rather than just the Cholesky factor of $\Sigma$), you can also condition on **structural shock paths**. For example, "assume the supply shock is zero for the next 4 periods" translates into direct constraints on elements of $\varepsilon$. + +Observable constraints and shock constraints can be combined in a single system. Observable constraints use the full MA representation (involving $\Phi$ and $P$), while shock constraints are simpler — they directly pin individual elements of $\varepsilon$. + +## Bayesian uncertainty + +The algorithm is applied independently to each posterior draw of $(B, \Sigma)$ or $(B, P)$. This means: + +- The unconditional forecast differs across draws (parameter uncertainty) +- The MA coefficients differ across draws +- The shock paths satisfying the constraints differ across draws + +The result is a full posterior distribution of conditional forecasts, from which you can compute medians, HDIs, and other summaries. The constrained variables will hit their targets exactly in every draw, but the unconstrained variables will show genuine posterior uncertainty about the conditional projection. diff --git a/docs/explanation/identification.md b/docs/explanation/identification.md index 5ff9b7f..2357730 100644 --- a/docs/explanation/identification.md +++ b/docs/explanation/identification.md @@ -37,3 +37,25 @@ scheme = SignRestriction( ``` Sign restrictions are weaker than Cholesky (they don't uniquely identify the model), but they require fewer assumptions. + +## Long-run restrictions (Blanchard-Quah) + +Cholesky and sign restrictions both constrain the **contemporaneous** (impact) response to structural shocks. Long-run restrictions take a different approach: they constrain the **cumulative** response as the horizon goes to infinity. + +The key object is the **long-run multiplier matrix**: + +$$C(1) = (I - A_1 - A_2 - \cdots - A_p)^{-1}$$ + +This matrix captures the total cumulative effect of a one-time shock. If the VAR is stationary, all shocks are transitory and $C(1)$ is finite. The long-run impact of structural shocks is $C(1) P$, where $P$ is the structural impact matrix. + +Blanchard & Quah (1989) proposed forcing $C(1) P$ to be **lower triangular**. This is achieved by applying the Cholesky decomposition not to the residual covariance $\Sigma$ (as in short-run Cholesky) but to the long-run covariance: + +$$C(1) \, \Sigma \, C(1)^\top = L \, L^\top$$ + +The structural impact matrix is then $P = C(1)^{-1} L$. + +The interpretation depends on the variable ordering: +- Shocks later in the ordering have **zero long-run cumulative effect** on variables earlier in the ordering +- The first shock is unrestricted in its long-run effects + +This is exactly the same Cholesky math, applied to a different matrix. The economics, however, are very different: you're restricting permanent effects rather than contemporaneous effects. In the classic Blanchard-Quah example, ordering output before prices means the second shock (interpreted as a demand shock) has no permanent effect on output — only supply shocks do. diff --git a/docs/explanation/minnesota-prior.md b/docs/explanation/minnesota-prior.md index 7e7348d..6e11369 100644 --- a/docs/explanation/minnesota-prior.md +++ b/docs/explanation/minnesota-prior.md @@ -1,6 +1,6 @@ # The Minnesota Prior -The **Minnesota prior** (Impulso, 1986) is the most widely used prior for Bayesian VARs. It encodes the belief that each variable follows a random walk, with coefficients on other variables' lags shrunk toward zero. +The **Minnesota prior** (Litterman, 1986) is the most widely used prior for Bayesian VARs. It encodes the belief that each variable follows a random walk, with coefficients on other variables' lags shrunk toward zero. ## Key hyperparameters @@ -27,3 +27,70 @@ spec = VAR(lags=4, prior="minnesota") prior = MinnesotaPrior(tightness=0.2, decay="geometric", cross_shrinkage=0.3) spec = VAR(lags=4, prior=prior) ``` + +## Conjugate estimation + +When the Minnesota prior is paired with a Normal-Inverse-Wishart (NIW) likelihood, the posterior belongs to the same family — this is called **conjugacy**. The practical consequence is dramatic: instead of running an iterative MCMC sampler, we can compute the posterior parameters in closed form and draw from it directly. + +The prior specifies: + +$$B \mid \Sigma \sim \mathcal{MN}(B_0, \Sigma, V_0), \qquad \Sigma \sim \mathcal{IW}(S_0, \nu_0)$$ + +where $\mathcal{MN}$ is the matrix normal distribution and $\mathcal{IW}$ is the inverse Wishart. Given data matrices $Y$ (observations) and $X$ (lagged regressors), the posterior parameters are: + +$$V_{\text{post}} = (V_0^{-1} + X^\top X)^{-1}$$ + +$$B_{\text{post}} = V_{\text{post}}(V_0^{-1} B_0 + X^\top Y)$$ + +$$\nu_{\text{post}} = \nu_0 + T$$ + +$$S_{\text{post}} = S_0 + Y^\top Y + B_0^\top V_0^{-1} B_0 - B_{\text{post}}^\top V_{\text{post}}^{-1} B_{\text{post}}$$ + +The posterior mean for $B$ is a precision-weighted average of the prior mean $B_0$ and the OLS estimate $(X^\top X)^{-1} X^\top Y$. When the prior is tight (small $V_0$), the posterior stays close to the random walk prior. When data is abundant (large $X^\top X$), the posterior approaches OLS. The posterior for $\Sigma$ is an inverse Wishart with updated scale and degrees of freedom. + +Sampling is straightforward: draw $\Sigma \sim \mathcal{IW}(S_{\text{post}}, \nu_{\text{post}})$, then $B \mid \Sigma \sim \mathcal{MN}(B_{\text{post}}, \Sigma, V_{\text{post}})$. Each draw is independent — no burn-in, no autocorrelation, no convergence diagnostics. `ConjugateVAR` implements this approach. + +## Data-driven prior selection + +The Minnesota prior's hyperparameters — `tightness`, `cross_shrinkage`, and `decay` — are often chosen by convention or trial-and-error. Giannone, Lenza & Primiceri (2015) proposed an empirical Bayes approach: choose hyperparameters by maximising the **marginal likelihood**. + +The marginal likelihood is the probability of the observed data $Y$ given hyperparameters $\lambda$, after integrating out all model parameters: + +$$p(Y \mid \lambda) = \int p(Y \mid B, \Sigma) \, p(B, \Sigma \mid \lambda) \, dB \, d\Sigma$$ + +Thanks to conjugacy, this integral has a closed-form solution involving determinants and multivariate gamma functions. Optimising $\log p(Y \mid \lambda)$ over $\lambda = (\text{tightness}, \text{cross\_shrinkage})$ using standard numerical optimisation (L-BFGS-B) is fast and reliable. + +The marginal likelihood naturally balances two forces: + +- **Fit**: a loose prior (high tightness) lets the model fit the data closely +- **Parsimony**: a loose prior also spreads probability mass over implausible parameter values, reducing the marginal likelihood + +The optimal hyperparameters sit at the sweet spot where the prior is just flexible enough to capture the data's patterns without wasting probability on noise. + +This approach is sometimes called **GLP** after the authors' initials. In Impulso, use `ConjugateVAR.optimize_prior()` or the shorthand `prior="minnesota_optimized"`. + +## Dummy observation priors + +An elegant trick from the classical VAR literature encodes prior beliefs not by modifying the prior distribution directly, but by **appending synthetic observations** to the dataset. These "dummy observations" push the posterior in the desired direction while preserving conjugacy. + +### Sum-of-coefficients prior + +The sum-of-coefficients prior (Doan, Litterman & Sims, 1984) encodes the belief that if all variables have been at their initial sample values forever, they should persist at those values. Formally, it adds observations implying: + +$$\sum_{j=1}^{p} A_j \approx I$$ + +where $A_j$ are the lag coefficient matrices. This discourages the model from predicting rapid mean-reversion when it isn't supported by the data. The hyperparameter $\mu$ controls the prior's strength: smaller $\mu$ imposes the belief more tightly. + +### Single-unit-root prior + +The single-unit-root prior (Sims, 1993) encodes the belief that each variable follows an independent random walk. It adds a single observation per variable that makes the model reluctant to introduce cointegrating relationships not strongly supported by the data. The hyperparameter $\delta$ controls its strength. + +### Practical usage + +In Impulso, dummy observations are added via `VARData.with_dummy_observations()`: + +```python +data_augmented = data.with_dummy_observations(n_lags=4, mu=1.0, delta=1.0) +``` + +The augmented data can then be passed to any estimation method — `ConjugateVAR`, `VAR`, or used with `optimize_prior()`. Because the dummy observations are simply extra rows in the data matrices, they preserve conjugacy and integrate seamlessly with marginal likelihood optimisation. diff --git a/docs/how-to/conditional-forecasting.md b/docs/how-to/conditional-forecasting.md new file mode 100644 index 0000000..18f0812 --- /dev/null +++ b/docs/how-to/conditional-forecasting.md @@ -0,0 +1,113 @@ +# Conditional Forecasting + +Standard forecasts let the model speak freely — given the historical data, where do the variables go next? But policy analysis often needs to answer a different question: **what happens if we assume a specific path for one variable?** + +For example: +- "What happens to GDP and inflation if the central bank raises the policy rate by 25bp at each of the next 4 meetings?" +- "What is the inflation outlook if oil prices stay at \$80/barrel for the next year?" + +Conditional forecasting answers these questions. It finds the **smallest set of shocks** consistent with the assumed path, then traces out the implications for all other variables. This implements the algorithm of Waggoner & Zha (1999). + +## Defining conditions + +A `ForecastCondition` specifies the variable, the periods (0-indexed forecast steps), and the target values: + +```python +from impulso import ForecastCondition + +# "The policy rate will be 5.25 at steps 0, 1, 2, and 3" +rate_path = ForecastCondition( + variable="rate", + periods=[0, 1, 2, 3], + values=[5.25, 5.50, 5.75, 6.00], +) +``` + +You can specify multiple conditions on different variables: + +```python +# Also condition on oil prices +oil_path = ForecastCondition( + variable="oil", + periods=[0, 1, 2, 3], + values=[80.0, 80.0, 80.0, 80.0], +) +``` + +!!! note "Periods are 0-indexed" + Period 0 is the first forecast step (one step ahead of the last observation). Period 3 is four steps ahead. + +## Reduced-form conditional forecasts + +On a `FittedVAR` (before identification), conditional forecasts use the Cholesky factor of the residual covariance to define the shock space: + +```python +from impulso import VAR, VARData, ForecastCondition + +data = VARData.from_df(df, endog=["gdp", "inflation", "rate"]) +fitted = VAR(lags=4, prior="minnesota").fit(data) + +# Condition on a rising rate path +rate_path = ForecastCondition( + variable="rate", + periods=[0, 1, 2, 3], + values=[5.25, 5.50, 5.75, 6.00], +) + +result = fitted.conditional_forecast(steps=8, conditions=[rate_path]) +``` + +The result is a `ConditionalForecastResult` — a subclass of `ForecastResult` that also stores the conditions. You can inspect it the same way: + +```python +# Posterior median conditional forecast +result.median() + +# HDI credible intervals +result.hdi(prob=0.89) + +# Plot +result.plot() +``` + +The constrained variable ("rate") will hit its target values exactly at the specified periods. The unconstrained variables ("gdp", "inflation") show the model's best prediction given those constraints — how the economy would evolve under the assumed rate path. + +!!! tip "Interpreting the results" + The conditional forecast answers: "What is the most likely path for all variables, given that `rate` follows the specified path?" The uncertainty bands for unconstrained variables reflect both parameter uncertainty (different posterior draws give different answers) and the uncertainty about which combination of shocks would produce the assumed path. + +## Structural conditional forecasts + +If you have an identified model, you can also condition on **structural shock paths**. This is useful for scenario analysis: "What happens if there are no supply shocks for the next 4 quarters?" + +```python +from impulso.identification import Cholesky + +identified = fitted.set_identification_strategy( + Cholesky(ordering=["gdp", "inflation", "rate"]) +) + +# Condition on zero supply shocks +no_supply_shocks = ForecastCondition( + variable="gdp", # shock named after the first variable in ordering + periods=[0, 1, 2, 3], + values=[0.0, 0.0, 0.0, 0.0], +) + +result = identified.conditional_forecast( + steps=8, + conditions=[rate_path], + shock_conditions=[no_supply_shocks], +) +``` + +Here, `conditions` constrain observable variable paths (as before), and `shock_conditions` constrain structural shock paths. You can use either or both. + +!!! warning "Shock naming" + Shock names correspond to the variable names in the identification scheme's ordering. With `Cholesky(ordering=["gdp", "inflation", "rate"])`, the shocks are named "gdp", "inflation", and "rate". + +## Degrees of freedom + +Each condition uses up one degree of freedom per constrained period. The total number of constraints cannot exceed `steps * n_variables` (the total number of shock values to be determined). In practice, keeping constraints well below this limit gives more stable results — the system becomes increasingly sensitive as you approach the maximum. + +!!! tip "Start simple" + Begin with conditions on one variable at a few periods. Add more constraints incrementally and check that results remain sensible. Over-constraining can produce large, implausible shocks. diff --git a/docs/how-to/conjugate-estimation.md b/docs/how-to/conjugate-estimation.md new file mode 100644 index 0000000..6505977 --- /dev/null +++ b/docs/how-to/conjugate-estimation.md @@ -0,0 +1,163 @@ +# Fast Estimation with ConjugateVAR + +Standard Bayesian VAR estimation uses Markov chain Monte Carlo (MCMC) — typically the NUTS sampler — to explore the posterior distribution of model parameters. This is flexible but slow: a moderate-sized model might take minutes to sample, and you need to worry about convergence diagnostics, burn-in, and chain autocorrelation. + +When your prior is Minnesota-type, none of this is necessary. The Minnesota prior combined with a Normal-Inverse-Wishart (NIW) likelihood gives a **conjugate** posterior — meaning the posterior has the same functional form as the prior. The posterior parameters can be computed in closed form, and draws are independent and identically distributed. No iteration, no burn-in, no convergence worries. + +`ConjugateVAR` implements this direct sampling approach. + +## Basic usage + +```python +from impulso import ConjugateVAR, VARData + +# Prepare your data as usual +data = VARData.from_df(df, endog=["gdp", "inflation", "rate"]) + +# Estimate — this returns a FittedVAR, just like VAR(...).fit() +model = ConjugateVAR(lags=4, prior="minnesota") +fitted = model.fit(data) +``` + +The resulting `FittedVAR` is identical to what you'd get from `VAR(...).fit(data, sampler)` — you can call `.forecast()`, `.set_identification_strategy()`, and everything else the same way. The only difference is how the posterior draws were obtained. + +!!! tip "When to use ConjugateVAR vs VAR" + Use `ConjugateVAR` when you're happy with a Minnesota-type prior and want speed. Use `VAR` with a NUTS sampler when you need custom priors, non-standard likelihoods, or stochastic volatility — things that break conjugacy. + +## Forecasting + +Since `ConjugateVAR.fit()` returns a standard `FittedVAR`, forecasting works exactly as before: + +```python +forecast = fitted.forecast(steps=8) + +# Posterior median forecast +forecast.median() + +# 89% highest density interval +forecast.hdi(prob=0.89) + +# Plot fan chart +forecast.plot() +``` + +The forecasts are probabilistic — each of the 2000 posterior draws (by default) produces a different forecast path. The median and HDI summarise this distribution. + +## Data-driven prior selection + +The Minnesota prior has hyperparameters — `tightness`, `cross_shrinkage`, and `decay` — that control how aggressively the posterior is pulled toward the prior mean. Choosing these by hand is common but somewhat arbitrary. + +**Giannone, Lenza & Primiceri (2015)** proposed a principled alternative: choose hyperparameters by maximising the **marginal likelihood** — the probability of the observed data given the hyperparameters, after integrating out all model parameters. This automatically balances fit and parsimony: too-tight priors underfit, too-loose priors overfit, and the marginal likelihood finds the sweet spot. + +Because the model is conjugate, the marginal likelihood has a closed-form expression — no additional sampling is needed. + +### One-step approach + +The simplest way is the `"minnesota_optimized"` shorthand, which optimises the prior and fits the model in one call: + +```python +fitted = ConjugateVAR(lags=4, prior="minnesota_optimized").fit(data) +``` + +This is equivalent to calling `optimize_prior()` followed by `fit()` with the optimised prior. + +### Two-step approach + +If you want to inspect or modify the optimised prior before fitting: + +```python +model = ConjugateVAR(lags=4) + +# Step 1: Find optimal hyperparameters +optimal_prior = model.optimize_prior(data) +print(optimal_prior) +# MinnesotaPrior(tightness=0.073, cross_shrinkage=0.42, decay='harmonic') + +# Step 2: Fit with the optimised prior +fitted = ConjugateVAR(lags=4, prior=optimal_prior).fit(data) +``` + +The optimiser adjusts `tightness` and `cross_shrinkage` while preserving the `decay` setting from the starting prior. + +!!! note "What the marginal likelihood measures" + The marginal likelihood answers: "how well does this combination of hyperparameters predict the observed data, averaging over all possible parameter values?" A higher value means the prior is better calibrated to the data. It naturally penalises both underfitting (prior too tight, can't match the data) and overfitting (prior too loose, wastes probability mass on implausible parameter values). + +### Comparing models + +You can also use the marginal likelihood to compare models with different lag orders: + +```python +for p in [1, 2, 3, 4]: + ml = ConjugateVAR(lags=p).marginal_likelihood(data) + print(f"Lags={p}: log ML = {ml:.1f}") +``` + +Higher log marginal likelihood indicates better fit after accounting for complexity. + +## Dummy observation priors + +Before fitting, you can augment your data with **dummy observations** that encode beliefs about persistence and unit roots. This is an elegant trick from the VAR literature (Doan, Litterman & Sims, 1984; Sims, 1993): rather than modifying the prior directly, you append synthetic data rows that push the posterior in the desired direction. + +### Sum-of-coefficients prior (mu) + +The sum-of-coefficients prior encodes the belief that **if all variables have been at their initial sample values forever, they should stay there**. This is a form of persistence belief — it discourages the model from predicting rapid mean-reversion that isn't supported by the data. + +The hyperparameter `mu` controls the strength: larger values mean a weaker prior (less influence on the posterior). + +```python +# Augment data with sum-of-coefficients dummy observations +data_augmented = data.with_dummy_observations(n_lags=4, mu=1.0) +``` + +### Single-unit-root prior (delta) + +The single-unit-root prior encodes the belief that **each variable follows a random walk** — unit root behaviour. This is related to cointegration beliefs and helps prevent spurious cointegrating relationships. + +The hyperparameter `delta` controls the strength: larger values mean a weaker prior. + +```python +# Both priors together +data_augmented = data.with_dummy_observations(n_lags=4, mu=1.0, delta=1.0) +``` + +### Combining with GLP optimisation + +Dummy observations work seamlessly with conjugate estimation and prior optimisation: + +```python +# Augment data, then optimise and fit +data_augmented = data.with_dummy_observations(n_lags=4, mu=1.0, delta=1.0) +fitted = ConjugateVAR(lags=4, prior="minnesota_optimized").fit(data_augmented) +``` + +!!! warning "Choosing mu and delta" + Start with `mu=1.0` and `delta=1.0` (moderate strength). Values below 0.5 impose strong beliefs; values above 5.0 have little effect. If you're unsure, the GLP marginal likelihood optimisation can help guide the choice — compare marginal likelihoods across different `mu`/`delta` combinations. + +## Complete workflow + +Putting it all together — dummy observations, optimised prior, estimation, and forecasting: + +```python +from impulso import ConjugateVAR, VARData +from impulso.identification import Cholesky + +# Prepare data +data = VARData.from_df(df, endog=["gdp", "inflation", "rate"]) + +# Augment with dummy observation priors +data = data.with_dummy_observations(n_lags=4, mu=1.0, delta=1.0) + +# Estimate with data-driven prior selection +fitted = ConjugateVAR(lags=4, prior="minnesota_optimized").fit(data) + +# Forecast +forecast = fitted.forecast(steps=8) +forecast.plot() + +# Structural analysis (works identically to VAR+NUTS path) +identified = fitted.set_identification_strategy( + Cholesky(ordering=["gdp", "inflation", "rate"]) +) +irf = identified.impulse_response(horizon=20) +irf.plot() +``` diff --git a/docs/how-to/long-run-restrictions.md b/docs/how-to/long-run-restrictions.md new file mode 100644 index 0000000..e56364d --- /dev/null +++ b/docs/how-to/long-run-restrictions.md @@ -0,0 +1,86 @@ +# Long-Run Restrictions + +Cholesky identification assumes a **contemporaneous** causal ordering: the first variable isn't affected by any other variable within the same period, the second is affected only by the first, and so on. This is a strong assumption that may not match your economic theory. + +Sometimes theory says nothing about contemporaneous effects but makes clear predictions about **long-run** effects. For example, in the Blanchard-Quah (1989) framework: + +- **Supply shocks** can have permanent effects on both output and prices +- **Demand shocks** have no permanent effect on output (only transitory) + +Long-run restrictions implement this idea. Instead of applying Cholesky to the contemporaneous impact matrix, we apply it to the **long-run cumulative impact matrix** — forcing some shocks to have zero permanent effect on certain variables. + +## Defining the scheme + +The ordering determines which shocks can have permanent effects: + +```python +from impulso.identification import LongRunRestriction + +scheme = LongRunRestriction(ordering=["output", "prices"]) +``` + +This means: +- The **first shock** (associated with "output") can have permanent effects on both output and prices +- The **second shock** (associated with "prices") has **no permanent effect on output** — only on prices + +The ordering encodes your identifying assumption about long-run neutrality. + +!!! tip "Reading the ordering" + Think of it as: "shocks later in the ordering cannot permanently affect variables earlier in the ordering." The last shock has the most restrictions; the first shock is unrestricted. + +## Applying to a fitted model + +The workflow is identical to Cholesky — only the economics differ: + +```python +from impulso import VAR, VARData +from impulso.identification import LongRunRestriction + +# Fit reduced-form VAR +data = VARData.from_df(df, endog=["output", "prices"]) +fitted = VAR(lags=4, prior="minnesota").fit(data) + +# Apply long-run identification +scheme = LongRunRestriction(ordering=["output", "prices"]) +identified = fitted.set_identification_strategy(scheme) +``` + +The resulting `IdentifiedVAR` supports all the same analysis methods — impulse responses, variance decomposition, and historical decomposition. + +## Impulse response analysis + +```python +irf = identified.impulse_response(horizon=40) +irf.plot() +``` + +When you examine the impulse responses, you should see the long-run restriction at work: the cumulative response of "output" to the second shock (the "prices" shock) converges to zero as the horizon increases. The first shock (the "output" shock) can have a permanent effect on both variables. + +!!! note "Convergence to zero" + The restriction forces the **cumulative** long-run effect to be zero, not the response at each individual horizon. The impulse response at horizon $h$ may be non-zero — it's only the sum across all horizons that vanishes. + +## Variance decomposition and historical decomposition + +These work identically to other identification schemes: + +```python +# What fraction of forecast error variance is due to each shock? +fevd = identified.fevd(horizon=40) +fevd.plot() + +# What drove the historical movements in each variable? +hd = identified.historical_decomposition() +hd.plot() +``` + +## When to use long-run restrictions + +| Situation | Recommended scheme | +|-----------|-------------------| +| Theory specifies contemporaneous ordering | `Cholesky` | +| Theory specifies long-run neutrality | `LongRunRestriction` | +| Theory specifies signs but not ordering | `SignRestriction` | +| Multiple schemes plausible | Try several and compare IRFs | + +!!! warning "Stationarity required" + Long-run restrictions require the VAR to be stationary (all eigenvalues of the companion matrix inside the unit circle). If the model has a unit root, the long-run multiplier matrix is undefined. If some posterior draws are near-non-stationary, results may be numerically unstable. diff --git a/docs/index.md b/docs/index.md index 95bd061..59b5c3a 100644 --- a/docs/index.md +++ b/docs/index.md @@ -2,41 +2,60 @@ **Bayesian Vector Autoregression in Python.** -```python -import pandas as pd -from impulso import VAR, VARData -from impulso.identification import Cholesky - -# Load data -df = pd.read_csv("macro_data.csv", index_col="date", parse_dates=True) -data = VARData.from_df(df, endog=["gdp", "inflation", "rate"]) - -# Estimate -fitted = VAR(lags="bic", prior="minnesota").fit(data) - -# Forecast -forecast = fitted.forecast(steps=8) -forecast.median() # point forecasts -forecast.hdi() # credible intervals - -# Structural analysis -identified = fitted.set_identification_strategy( - Cholesky(ordering=["gdp", "inflation", "rate"]) -) -irf = identified.impulse_response(horizon=20) -irf.plot() -``` +=== "Fast (Conjugate)" + + ```python + import pandas as pd + from impulso import ConjugateVAR, VARData + + # Load data + df = pd.read_csv("macro_data.csv", index_col="date", parse_dates=True) + data = VARData.from_df(df, endog=["gdp", "inflation", "rate"]) + + # Estimate with data-driven prior selection (no MCMC needed) + fitted = ConjugateVAR(lags=4, prior="minnesota_optimized").fit(data) + + # Forecast + forecast = fitted.forecast(steps=8) + forecast.median() # point forecasts + forecast.hdi() # credible intervals + ``` + +=== "Flexible (NUTS)" + + ```python + import pandas as pd + from impulso import VAR, VARData + from impulso.identification import Cholesky + + # Load data + df = pd.read_csv("macro_data.csv", index_col="date", parse_dates=True) + data = VARData.from_df(df, endog=["gdp", "inflation", "rate"]) + + # Estimate with NUTS (supports arbitrary priors) + fitted = VAR(lags="bic", prior="minnesota").fit(data) + + # Structural analysis + identified = fitted.set_identification_strategy( + Cholesky(ordering=["gdp", "inflation", "rate"]) + ) + irf = identified.impulse_response(horizon=20) + irf.plot() + ``` ## Features - **Validated data containers** — `VARData` catches shape mismatches, missing values, and type errors at construction time - **Immutable pipeline** — `VARData` -> `VAR` -> `FittedVAR` -> `IdentifiedVAR`, each stage frozen after creation - **Economist-friendly API** — think in variables and lags, not tensors and MCMC chains +- **Two estimation paths** — conjugate NIW sampling (instant, no MCMC) or full NUTS via PyMC (flexible, supports custom priors) +- **Data-driven prior selection** — automatic Minnesota hyperparameter optimisation via marginal likelihood (Giannone, Lenza & Primiceri, 2015) +- **Dummy observation priors** — encode persistence and unit root beliefs via sum-of-coefficients and single-unit-root priors - **Minnesota prior** — smart defaults with tunable hyperparameters for shrinkage - **Automatic lag selection** — AIC, BIC, and Hannan-Quinn criteria -- **PyMC backend** — full Bayesian estimation with NUTS sampling - **Probabilistic forecasts** — posterior median, HDI credible intervals, tidy DataFrames -- **Structural identification** — Cholesky and sign restriction schemes +- **Conditional forecasting** — constrain future variable paths or structural shock paths for policy analysis +- **Structural identification** — Cholesky, sign restriction, and Blanchard-Quah long-run schemes - **Built-in plotting** — IRF, FEVD, forecast, and historical decomposition plots ## Installation @@ -47,7 +66,7 @@ pip install impulso ## Learn more -- [Quickstart tutorial](tutorials/quickstart.ipynb) — fit your first Bayesian VAR -- [Forecasting tutorial](tutorials/forecasting.ipynb) — produce probabilistic forecasts -- [Structural analysis tutorial](tutorials/structural-analysis.ipynb) — impulse responses and variance decompositions +- [Fast estimation with ConjugateVAR](how-to/conjugate-estimation.md) — conjugate sampling, data-driven priors, dummy observations +- [Conditional forecasting](how-to/conditional-forecasting.md) — constrain future paths for policy analysis +- [Long-run restrictions](how-to/long-run-restrictions.md) — Blanchard-Quah structural identification - [API Reference](reference/index.md) — complete module documentation diff --git a/docs/reference/conditions.md b/docs/reference/conditions.md new file mode 100644 index 0000000..bc68fb0 --- /dev/null +++ b/docs/reference/conditions.md @@ -0,0 +1,3 @@ +# Forecast Conditions + +::: impulso.conditions diff --git a/docs/reference/conjugate.md b/docs/reference/conjugate.md new file mode 100644 index 0000000..11358e9 --- /dev/null +++ b/docs/reference/conjugate.md @@ -0,0 +1,3 @@ +# ConjugateVAR + +::: impulso.conjugate diff --git a/mkdocs.yml b/mkdocs.yml index 536707e..099169d 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -11,8 +11,6 @@ nav: - Home: index.md - Tutorials: - tutorials/index.md - # - tutorials/quickstart.ipynb - # - tutorials/forecasting.ipynb - tutorials/structural-analysis.ipynb - How-To Guides: - how-to/index.md @@ -20,15 +18,21 @@ nav: - how-to/custom-priors.md - how-to/lag-selection.md - how-to/sign-restrictions.md + - how-to/conjugate-estimation.md + - how-to/long-run-restrictions.md + - how-to/conditional-forecasting.md - Explanation: - explanation/index.md - explanation/bayesian-var.md - explanation/minnesota-prior.md - explanation/identification.md + - explanation/conditional-forecasting.md - Reference: - reference/index.md - reference/data.md - reference/spec.md + - reference/conjugate.md + - reference/conditions.md - reference/priors.md - reference/samplers.md - reference/fitted.md diff --git a/src/impulso/__init__.py b/src/impulso/__init__.py index 790f00b..31f5a57 100644 --- a/src/impulso/__init__.py +++ b/src/impulso/__init__.py @@ -7,14 +7,18 @@ __all__ = [ "VAR", "Cholesky", + "ConditionalForecastResult", + "ConjugateVAR", "FEVDResult", "FittedVAR", + "ForecastCondition", "ForecastResult", "HDIResult", "HistoricalDecompositionResult", "IRFResult", "IdentifiedVAR", "LagOrderResult", + "LongRunRestriction", "MinnesotaPrior", "NUTSSampler", "SignRestriction", @@ -27,9 +31,13 @@ def __getattr__(name: str): """Lazy imports for types not needed at import time.""" _lazy_imports = { + "ConditionalForecastResult": "impulso.results", + "ConjugateVAR": "impulso.conjugate", "FittedVAR": "impulso.fitted", + "ForecastCondition": "impulso.conditions", "IdentifiedVAR": "impulso.identified", "Cholesky": "impulso.identification", + "LongRunRestriction": "impulso.identification", "SignRestriction": "impulso.identification", "MinnesotaPrior": "impulso.priors", "NUTSSampler": "impulso.samplers", diff --git a/src/impulso/conditions.py b/src/impulso/conditions.py new file mode 100644 index 0000000..06f438d --- /dev/null +++ b/src/impulso/conditions.py @@ -0,0 +1,33 @@ +"""Forecast condition definitions for conditional forecasting.""" + +from typing import Literal, Self + +from pydantic import model_validator + +from impulso._base import ImpulsoModel + + +class ForecastCondition(ImpulsoModel): + """A constraint on a variable's future path for conditional forecasting. + + Attributes: + variable: Name of the variable to constrain. + periods: Forecast steps to constrain (0-indexed). + values: Target values at those periods. + constraint_type: Type of constraint. Only 'hard' is currently supported. + """ + + variable: str + periods: list[int] + values: list[float] + constraint_type: Literal["hard"] = "hard" + + @model_validator(mode="after") + def _validate_periods_values_match(self) -> Self: + if len(self.periods) != len(self.values): + raise ValueError(f"periods length ({len(self.periods)}) must equal values length ({len(self.values)})") + if len(self.periods) == 0: + raise ValueError("periods must be non-empty") + if any(p < 0 for p in self.periods): + raise ValueError("All periods must be non-negative") + return self diff --git a/src/impulso/conjugate.py b/src/impulso/conjugate.py new file mode 100644 index 0000000..ea59ab8 --- /dev/null +++ b/src/impulso/conjugate.py @@ -0,0 +1,348 @@ +"""ConjugateVAR — direct Normal-Inverse-Wishart posterior sampling.""" + +from typing import TYPE_CHECKING, Literal, Self + +import numpy as np +from pydantic import Field, model_validator + +from impulso._base import ImpulsoBaseModel +from impulso.data import VARData +from impulso.priors import MinnesotaPrior + +if TYPE_CHECKING: + from impulso.fitted import FittedVAR + + +class ConjugateVAR(ImpulsoBaseModel): + """Bayesian VAR with conjugate Normal-Inverse-Wishart estimation. + + Produces iid posterior draws via direct sampling — no MCMC iteration, + no burn-in, no autocorrelation. Orders of magnitude faster than NUTS + for models with Minnesota-type priors. + + Attributes: + lags: Fixed lag order or selection criterion. + max_lags: Upper bound for automatic lag selection. + prior: Minnesota prior instance or string shorthand. + draws: Number of posterior draws. + random_seed: Seed for reproducibility. + """ + + lags: int | Literal["aic", "bic", "hq"] = Field(...) + max_lags: int | None = None + prior: Literal["minnesota", "minnesota_optimized"] | MinnesotaPrior = "minnesota" + draws: int = Field(2000, ge=1) + random_seed: int | None = None + + @model_validator(mode="after") + def _validate_spec(self) -> Self: + if self.max_lags is not None and isinstance(self.lags, int): + raise ValueError("max_lags is only valid when lags is a selection criterion") + if isinstance(self.lags, int) and self.lags < 1: + raise ValueError(f"lags must be >= 1, got {self.lags}") + return self + + @property + def resolved_prior(self) -> MinnesotaPrior: + """Resolve string shorthand to a MinnesotaPrior instance.""" + if isinstance(self.prior, str): + return MinnesotaPrior() + return self.prior + + @staticmethod + def _build_niw_params( + data: VARData, + n_lags: int, + prior: MinnesotaPrior, + ) -> dict[str, np.ndarray]: + """Build data matrices and compute NIW prior/posterior parameters. + + Args: + data: VARData instance. + n_lags: Number of lags. + prior: MinnesotaPrior instance. + + Returns: + Dictionary with keys: Y, X, B_prior, V_prior, V_prior_inv, + V_posterior, B_posterior, S_prior, S_posterior, nu_prior, nu_posterior. + """ + n_vars = data.endog.shape[1] + prior_params = prior.build_priors(n_vars=n_vars, n_lags=n_lags) + + # Build data matrices: Y = (T-p, n), X = (T-p, n*p + 1) with intercept + y = data.endog + Y = y[n_lags:] + X_parts = [np.ones((Y.shape[0], 1))] + for lag in range(1, n_lags + 1): + X_parts.append(y[n_lags - lag : -lag]) + X = np.hstack(X_parts) + + T_eff = Y.shape[0] + n_coeffs = X.shape[1] + + # Convert Minnesota prior to NIW parameters + # Prior mean: [intercept_prior | B_mu] + B_prior = np.zeros((n_coeffs, n_vars)) + B_prior[1:, :] = prior_params["B_mu"].T + + # Prior covariance: diagonal from B_sigma + # Intercept gets a wide prior (variance=1.0, matching PyMC path) + intercept_var = 1.0 + lag_var = np.mean(prior_params["B_sigma"] ** 2, axis=0) + prior_var_diag = np.concatenate([[intercept_var], lag_var]) + V_prior = np.diag(prior_var_diag) + V_prior_inv = np.diag(1.0 / prior_var_diag) + + # OLS estimates for scale matrix initialisation + B_ols = np.linalg.lstsq(X, Y, rcond=None)[0] + resid_ols = Y - X @ B_ols + sigma_ols = (resid_ols.T @ resid_ols) / T_eff + + # NIW prior hyperparameters + nu_prior = n_vars + 2 # minimally informative + S_prior = sigma_ols * (nu_prior - n_vars - 1) # centres IW mode at sigma_ols + + # Posterior parameters + V_posterior = np.linalg.inv(V_prior_inv + X.T @ X) + B_posterior = V_posterior @ (V_prior_inv @ B_prior + X.T @ Y) + nu_posterior = nu_prior + T_eff + S_posterior = ( + S_prior + + Y.T @ Y + + B_prior.T @ V_prior_inv @ B_prior + - B_posterior.T @ np.linalg.inv(V_posterior) @ B_posterior + ) + S_posterior = (S_posterior + S_posterior.T) / 2 + + return { + "Y": Y, + "X": X, + "B_prior": B_prior, + "V_prior": V_prior, + "V_prior_inv": V_prior_inv, + "V_posterior": V_posterior, + "B_posterior": B_posterior, + "S_prior": S_prior, + "S_posterior": S_posterior, + "nu_prior": nu_prior, + "nu_posterior": nu_posterior, + } + + def fit(self, data: VARData) -> "FittedVAR": + """Estimate the Bayesian VAR via conjugate NIW posterior sampling. + + Args: + data: VARData instance. + + Returns: + FittedVAR with iid posterior draws. + """ + import arviz as az + import xarray as xr + from scipy.stats import invwishart + + from impulso._lag_selection import select_lag_order + from impulso.fitted import FittedVAR + + # Resolve lags + if isinstance(self.lags, str): + max_lags = self.max_lags or 12 + ic = select_lag_order(data, max_lags=max_lags) + n_lags = getattr(ic, self.lags) + else: + n_lags = self.lags + + n_vars = data.endog.shape[1] + + # Resolve prior (optimize if requested) + if isinstance(self.prior, str) and self.prior == "minnesota_optimized": + prior = self._optimize_prior_internal(data, n_lags) + else: + prior = self.resolved_prior + + params = self._build_niw_params(data, n_lags, prior) + V_posterior = params["V_posterior"] + B_posterior = params["B_posterior"] + S_posterior = params["S_posterior"] + nu_posterior = params["nu_posterior"] + n_coeffs = B_posterior.shape[0] # 1 + n_vars * n_lags + + # Direct sampling + rng = np.random.default_rng(self.random_seed) + + B_draws = np.zeros((self.draws, n_coeffs, n_vars)) + Sigma_draws = np.zeros((self.draws, n_vars, n_vars)) + + chol_V_posterior = np.linalg.cholesky(V_posterior) + + for i in range(self.draws): + # Draw Sigma ~ IW(S_posterior, nu_posterior) + Sigma_draw = invwishart.rvs(df=nu_posterior, scale=S_posterior, random_state=rng) + Sigma_draws[i] = Sigma_draw + + # Draw B | Sigma ~ MN(B_posterior, Sigma, V_posterior) + # vec(B) ~ N(vec(B_posterior), Sigma kron V_posterior) + chol_Sigma = np.linalg.cholesky(Sigma_draw) + Z = rng.standard_normal((n_coeffs, n_vars)) + B_draw = B_posterior + chol_V_posterior @ Z @ chol_Sigma.T + B_draws[i] = B_draw + + # Separate intercept and lag coefficients + intercept_arr = B_draws[:, 0, :] # (draws, n_vars) + B_lag_arr = B_draws[:, 1:, :] # (draws, n*p, n_vars) + # Transpose to match PyMC convention: B is (n_vars, n_vars*n_lags) + B_lag_arr = np.swapaxes(B_lag_arr, -2, -1) # (draws, n_vars, n*p) + + # Add chain dimension (chains=1 for conjugate) + intercept_arr = intercept_arr[np.newaxis, :] # (1, draws, n_vars) + B_lag_arr = B_lag_arr[np.newaxis, :] # (1, draws, n_vars, n*p) + Sigma_draws = Sigma_draws[np.newaxis, :] # (1, draws, n_vars, n_vars) + + # Package as InferenceData + posterior = xr.Dataset({ + "B": xr.DataArray(B_lag_arr, dims=["chain", "draw", "equations", "coefficients"]), + "intercept": xr.DataArray(intercept_arr, dims=["chain", "draw", "equations"]), + "Sigma": xr.DataArray(Sigma_draws, dims=["chain", "draw", "var1", "var2"]), + }) + idata = az.InferenceData(posterior=posterior) + + return FittedVAR.model_construct( + idata=idata, + n_lags=n_lags, + data=data, + var_names=data.endog_names, + ) + + def optimize_prior( + self, + data: VARData, + ) -> MinnesotaPrior: + """Find Minnesota hyperparameters maximising the marginal likelihood. + + Implements Giannone, Lenza & Primiceri (2015) data-driven prior + selection via closed-form marginal likelihood optimisation. + + Args: + data: VARData instance (may include dummy observations). + + Returns: + MinnesotaPrior with optimal tightness and cross_shrinkage. + """ + from impulso._lag_selection import select_lag_order + + # Resolve lags + if isinstance(self.lags, str): + max_lags = self.max_lags or 12 + ic = select_lag_order(data, max_lags=max_lags) + n_lags = getattr(ic, self.lags) + else: + n_lags = self.lags + + return self._optimize_prior_internal(data, n_lags) + + def _optimize_prior_internal( + self, + data: VARData, + n_lags: int, + ) -> MinnesotaPrior: + """Internal implementation of prior optimisation.""" + from scipy.optimize import minimize + + current_prior = self.resolved_prior + + def neg_log_marginal_likelihood(params: np.ndarray) -> float: + tightness = params[0] + cross_shrinkage = params[1] + prior = MinnesotaPrior( + tightness=tightness, + cross_shrinkage=cross_shrinkage, + decay=current_prior.decay, + ) + return -self._log_marginal_likelihood(data, n_lags, prior) + + x0 = np.array([current_prior.tightness, current_prior.cross_shrinkage]) + bounds = [(0.001, 10.0), (0.01, 1.0)] + + result = minimize( + neg_log_marginal_likelihood, + x0=x0, + method="L-BFGS-B", + bounds=bounds, + ) + + return MinnesotaPrior( + tightness=float(result.x[0]), + cross_shrinkage=float(result.x[1]), + decay=current_prior.decay, + ) + + def _log_marginal_likelihood( + self, + data: VARData, + n_lags: int, + prior: MinnesotaPrior, + ) -> float: + """Compute log marginal likelihood p(Y|lambda) for NIW conjugate model. + + Uses the closed-form NIW marginal likelihood from Kadiyala & Karlsson (1997). + + Args: + data: VARData instance. + n_lags: Number of lags. + prior: MinnesotaPrior with specific hyperparameters. + + Returns: + Log marginal likelihood (scalar). + """ + from scipy.special import gammaln + + n_vars = data.endog.shape[1] + params = self._build_niw_params(data, n_lags, prior) + + T_eff = params["Y"].shape[0] + V_prior = params["V_prior"] + V_posterior = params["V_posterior"] + S_prior = params["S_prior"] + S_posterior = params["S_posterior"] + nu_prior = params["nu_prior"] + nu_posterior = params["nu_posterior"] + + # Log marginal likelihood formula + log_ml = 0.0 + log_ml -= (T_eff * n_vars / 2) * np.log(np.pi) + + # Log-determinant terms + _, logdet_V_prior = np.linalg.slogdet(V_prior) + _, logdet_V_posterior = np.linalg.slogdet(V_posterior) + log_ml += 0.5 * (logdet_V_posterior - logdet_V_prior) * n_vars + + _, logdet_S_prior = np.linalg.slogdet(S_prior) + _, logdet_S_posterior = np.linalg.slogdet(S_posterior) + log_ml += (nu_prior / 2) * logdet_S_prior + log_ml -= (nu_posterior / 2) * logdet_S_posterior + + # Multivariate gamma function terms + for j in range(n_vars): + log_ml += gammaln((nu_posterior - j) / 2) - gammaln((nu_prior - j) / 2) + + return log_ml + + def marginal_likelihood(self, data: VARData) -> float: + """Compute log marginal likelihood for the current prior. + + Args: + data: VARData instance. + + Returns: + Log marginal likelihood (scalar). + """ + from impulso._lag_selection import select_lag_order + + if isinstance(self.lags, str): + max_lags = self.max_lags or 12 + ic = select_lag_order(data, max_lags=max_lags) + n_lags = getattr(ic, self.lags) + else: + n_lags = self.lags + + return self._log_marginal_likelihood(data, n_lags, self.resolved_prior) diff --git a/src/impulso/data.py b/src/impulso/data.py index cb7fe54..62d229e 100644 --- a/src/impulso/data.py +++ b/src/impulso/data.py @@ -99,3 +99,67 @@ def from_df( exog_names=exog, index=df.index, ) + + def with_dummy_observations( + self, + n_lags: int, + mu: float | None = None, + delta: float | None = None, + ) -> "VARData": + """Return new VARData with dummy observations appended. + + Dummy observations encode beliefs about unit roots and persistence, + following Doan, Litterman & Sims (1984) and Sims (1993). + + Args: + n_lags: Number of VAR lags. Validated here for API consistency; + the dummy values themselves do not depend on n_lags, but + downstream fitting requires it to be valid. + mu: Sum-of-coefficients hyperparameter. Larger = weaker prior. + delta: Single-unit-root hyperparameter. Larger = weaker prior. + + Returns: + New VARData with dummy observations appended to endog. + """ + if mu is None and delta is None: + raise ValueError("At least one of mu or delta must be provided") + if mu is not None and mu <= 0: + raise ValueError(f"mu must be strictly positive, got {mu}") + if delta is not None and delta <= 0: + raise ValueError(f"delta must be strictly positive, got {delta}") + if n_lags < 1: + raise ValueError(f"n_lags must be >= 1, got {n_lags}") + + n_vars = self.endog.shape[1] + y_bar = self.endog.mean(axis=0) + dummy_rows = [] + + if mu is not None: + soc = np.zeros((n_vars, n_vars)) + np.fill_diagonal(soc, y_bar / mu) + dummy_rows.append(soc) + + if delta is not None: + sur = (y_bar / delta).reshape(1, n_vars) + dummy_rows.append(sur) + + dummies = np.vstack(dummy_rows) + new_endog = np.vstack([self.endog, dummies]) + + freq = self.index.freq or pd.tseries.frequencies.to_offset(pd.infer_freq(self.index)) + n_dummy = dummies.shape[0] + extra_index = pd.date_range(start=self.index[-1] + freq, periods=n_dummy, freq=freq) + new_index = self.index.append(extra_index) + + new_exog = None + if self.exog is not None: + exog_padding = np.zeros((n_dummy, self.exog.shape[1])) + new_exog = np.vstack([self.exog, exog_padding]) + + return VARData( + endog=new_endog, + endog_names=self.endog_names, + exog=new_exog, + exog_names=self.exog_names, + index=new_index, + ) diff --git a/src/impulso/fitted.py b/src/impulso/fitted.py index 9b16141..8dc6a88 100644 --- a/src/impulso/fitted.py +++ b/src/impulso/fitted.py @@ -11,8 +11,9 @@ from impulso.protocols import IdentificationScheme if TYPE_CHECKING: + from impulso.conditions import ForecastCondition from impulso.identified import IdentifiedVAR - from impulso.results import ForecastResult + from impulso.results import ConditionalForecastResult, ForecastResult class FittedVAR(ImpulsoBaseModel): @@ -111,6 +112,171 @@ def forecast( return ForecastResult(idata=idata, steps=steps, var_names=self.var_names) + @staticmethod + def _ma_coefficients_single(B: np.ndarray, n_vars: int, n_lags: int, steps: int) -> list[np.ndarray]: + """Compute MA coefficient recursion for a single draw. + + Args: + B: Coefficient matrix (n_vars, n_vars*n_lags). + n_vars: Number of endogenous variables. + n_lags: Number of lags. + steps: Number of forecast steps. + + Returns: + List of MA coefficient matrices [Phi_0, ..., Phi_{steps-1}]. + """ + A_matrices = [B[:, j * n_vars : (j + 1) * n_vars] for j in range(n_lags)] + ma_coefficients: list[np.ndarray] = [np.eye(n_vars)] + for h in range(1, steps): + phi_h = np.zeros((n_vars, n_vars)) + for j in range(min(h, n_lags)): + phi_h += A_matrices[j] @ ma_coefficients[h - j - 1] + ma_coefficients.append(phi_h) + return ma_coefficients + + @staticmethod + def _build_constraint_system( + conditions: "list[ForecastCondition]", + var_names: list[str], + ma_coefficients: list[np.ndarray], + impact_matrix: np.ndarray, + unconditional: np.ndarray, + steps: int, + n_vars: int, + ) -> tuple[np.ndarray, np.ndarray]: + """Build the linear constraint system R @ shocks = target. + + Args: + conditions: List of ForecastCondition instances. + var_names: Variable names. + ma_coefficients: MA coefficient matrices. + impact_matrix: Cholesky factor of Sigma or structural matrix P. + unconditional: Unconditional forecast array (steps, n_vars). + steps: Number of forecast steps. + n_vars: Number of variables. + + Returns: + Tuple of (R, target) arrays for the constraint system. + """ + constraint_rows = [] + constraint_targets = [] + for cond in conditions: + var_idx = var_names.index(cond.variable) + for period, value in zip(cond.periods, cond.values, strict=True): + row = np.zeros(steps * n_vars) + for s in range(period + 1): + response = ma_coefficients[period - s] @ impact_matrix + row[s * n_vars : (s + 1) * n_vars] = response[var_idx, :] + constraint_rows.append(row) + constraint_targets.append(value - unconditional[period, var_idx]) + return np.array(constraint_rows), np.array(constraint_targets) + + def conditional_forecast( + self, + steps: int, + conditions: "list[ForecastCondition]", + exog_future: np.ndarray | None = None, + ) -> "ConditionalForecastResult": + """Produce conditional forecasts subject to constraints on future paths. + + Implements the Waggoner & Zha (1999) algorithm for hard constraints. + Computes unconditional forecasts, then solves for the shock paths + that satisfy the constraints. + + Args: + steps: Number of forecast steps. + conditions: List of ForecastCondition instances specifying constraints. + exog_future: Future exogenous values if model has exog. + + Returns: + ConditionalForecastResult with constrained posterior forecast draws. + """ + import xarray as xr + + from impulso.results import ConditionalForecastResult + + self._validate_conditions(conditions, steps) + + B_draws = self.coefficients # (C, D, n, n*p) + intercept_draws = self.intercepts # (C, D, n) + sigma_draws = self.sigma # (C, D, n, n) + n_chains, n_draws, n_vars, _ = B_draws.shape + + y_hist = self.data.endog[-self.n_lags :] + forecasts = np.zeros((n_chains, n_draws, steps, n_vars)) + + for c in range(n_chains): + for d in range(n_draws): + B = B_draws[c, d] + intercept = intercept_draws[c, d] + chol_sigma = np.linalg.cholesky(sigma_draws[c, d]) + + ma_coefficients = self._ma_coefficients_single(B, n_vars, self.n_lags, steps) + unconditional = self._unconditional_forecast_single( + B, intercept, y_hist, steps, self.n_lags, n_vars, c, d, exog_future + ) + + R, target = self._build_constraint_system( + conditions, self.var_names, ma_coefficients, chol_sigma, unconditional, steps, n_vars + ) + + shocks, _, _, _ = np.linalg.lstsq(R, target, rcond=None) + shocks = shocks.reshape(steps, n_vars) + + conditional = unconditional.copy() + for h in range(steps): + for s in range(h + 1): + conditional[h] += ma_coefficients[h - s] @ chol_sigma @ shocks[s] + + forecasts[c, d] = conditional + + forecast_da = xr.DataArray( + forecasts, + dims=["chain", "draw", "step", "variable"], + coords={"variable": self.var_names}, + name="forecast", + ) + idata = az.InferenceData(posterior_predictive=xr.Dataset({"forecast": forecast_da})) + return ConditionalForecastResult(idata=idata, steps=steps, var_names=self.var_names, conditions=conditions) + + def _validate_conditions(self, conditions: "list[ForecastCondition]", steps: int) -> None: + """Validate forecast conditions against model variables and step range.""" + for cond in conditions: + if cond.variable not in self.var_names: + raise ValueError(f"Condition variable '{cond.variable}' not in var_names {self.var_names}") + for p in cond.periods: + if p < 0 or p >= steps: + raise ValueError(f"Condition period {p} out of range for {steps} forecast steps") + + def _unconditional_forecast_single( + self, + B: np.ndarray, + intercept: np.ndarray, + y_hist: np.ndarray, + steps: int, + n_lags: int, + n_vars: int, + chain_idx: int, + draw_idx: int, + exog_future: np.ndarray | None, + ) -> np.ndarray: + """Compute unconditional forecast for a single posterior draw. + + Returns: + Array of shape (steps, n_vars). + """ + y_buffer = y_hist.copy() + unconditional = np.zeros((steps, n_vars)) + for h in range(steps): + x_lag = np.concatenate([y_buffer[-(lag + 1)] for lag in range(n_lags)]) + y_new = intercept + B @ x_lag + if self.has_exog and exog_future is not None: + B_exog = self.idata.posterior["B_exog"].values[chain_idx, draw_idx] + y_new = y_new + B_exog @ exog_future[h] + unconditional[h] = y_new + y_buffer = np.vstack([y_buffer[1:], y_new.reshape(1, -1)]) + return unconditional + def set_identification_strategy(self, scheme: IdentificationScheme) -> "IdentifiedVAR": """Apply a structural identification scheme. diff --git a/src/impulso/identification.py b/src/impulso/identification.py index a6675cb..f03312b 100644 --- a/src/impulso/identification.py +++ b/src/impulso/identification.py @@ -192,3 +192,85 @@ def _check_restrictions(self, candidate: np.ndarray, var_names: list[str], shock if sign == "-" and val > 0: return False return True + + +class LongRunRestriction(ImpulsoModel): + """Blanchard-Quah long-run identification scheme. + + Identifies structural shocks by their long-run cumulative effects. + The long-run impact matrix is forced to be lower triangular via + Cholesky decomposition, so the first shock has no permanent effect + on the second variable, etc. + + Attributes: + ordering: Variable ordering (determines which shocks have + permanent effects on which variables). + """ + + ordering: list[str] + + def identify(self, idata: az.InferenceData, var_names: list[str]) -> az.InferenceData: + """Apply Blanchard-Quah long-run identification. + + Args: + idata: InferenceData with 'B' and 'Sigma' in posterior. + var_names: Variable names from the VAR model. + + Returns: + InferenceData with 'structural_shock_matrix' added to posterior. + """ + if "B" not in idata.posterior: + raise ValueError("LongRunRestriction requires 'B' (VAR coefficients) in idata.posterior") + + unknown = set(self.ordering) - set(var_names) + if unknown: + raise ValueError(f"ordering contains unknown variables: {unknown}") + if len(self.ordering) != len(var_names): + raise ValueError(f"ordering must contain exactly {len(var_names)} variables, got {len(self.ordering)}") + + B_draws = idata.posterior["B"].values # (C, D, n, n*p) + sigma_draws = idata.posterior["Sigma"].values # (C, D, n, n) + n_chains, n_draws, n_vars, n_total_coeffs = B_draws.shape + n_lags = n_total_coeffs // n_vars + + # Compute permutation for reordering + perm = [var_names.index(v) for v in self.ordering] + + P = np.zeros((n_chains, n_draws, n_vars, n_vars)) + + for c in range(n_chains): + for d in range(n_draws): + B = B_draws[c, d] # (n, n*p) + Sigma = sigma_draws[c, d] # (n, n) + + # Sum of lag coefficient matrices: A_1 + A_2 + ... + A_p + lag_coefficient_sum = np.zeros((n_vars, n_vars)) + for j in range(n_lags): + lag_coefficient_sum += B[:, j * n_vars : (j + 1) * n_vars] + + # Long-run multiplier: (I - A_1 - ... - A_p)^{-1} + long_run_multiplier = np.linalg.inv(np.eye(n_vars) - lag_coefficient_sum) + + # Reorder for requested ordering + long_run_multiplier_ordered = long_run_multiplier[np.ix_(perm, perm)] + Sigma_ordered = Sigma[np.ix_(perm, perm)] + + # Long-run covariance + long_run_covariance = long_run_multiplier_ordered @ Sigma_ordered @ long_run_multiplier_ordered.T + + # Cholesky of long-run covariance + long_run_cholesky = np.linalg.cholesky(long_run_covariance) + + # Structural impact matrix: P = C(1)^{-1} @ L + structural_impact_matrix = np.linalg.inv(long_run_multiplier_ordered) @ long_run_cholesky + + P[c, d] = structural_impact_matrix + + P_da = xr.DataArray( + P, + dims=["chain", "draw", "response", "shock"], + coords={"response": self.ordering, "shock": self.ordering}, + ) + + new_posterior = idata.posterior.assign(structural_shock_matrix=P_da) + return az.InferenceData(posterior=new_posterior) diff --git a/src/impulso/identified.py b/src/impulso/identified.py index b94dc82..911937a 100644 --- a/src/impulso/identified.py +++ b/src/impulso/identified.py @@ -1,5 +1,7 @@ """IdentifiedVAR — structural VAR with identified shocks.""" +from typing import TYPE_CHECKING + import arviz as az import numpy as np import pandas as pd @@ -10,6 +12,10 @@ from impulso.data import VARData from impulso.results import FEVDResult, HistoricalDecompositionResult, IRFResult +if TYPE_CHECKING: + from impulso.conditions import ForecastCondition + from impulso.results import ConditionalForecastResult + class IdentifiedVAR(ImpulsoBaseModel): """Immutable structural VAR with identified shocks. @@ -179,3 +185,194 @@ def historical_decomposition( ) idata = az.InferenceData(posterior_predictive=xr.Dataset({"hd": hd_da})) return HistoricalDecompositionResult(idata=idata, var_names=self.var_names) + + def conditional_forecast( + self, + steps: int, + conditions: "list[ForecastCondition]", + shock_conditions: "list[ForecastCondition] | None" = None, + exog_future: np.ndarray | None = None, + ) -> "ConditionalForecastResult": + """Produce structural conditional forecasts. + + Extends reduced-form conditional forecasting by allowing constraints + on structural shock paths in addition to observable variable paths. + + Args: + steps: Number of forecast steps. + conditions: List of ForecastCondition instances for observables. + shock_conditions: Optional list of ForecastCondition instances for + structural shocks. + exog_future: Future exogenous values if model has exog. + + Returns: + ConditionalForecastResult with constrained forecast draws. + """ + from impulso.fitted import FittedVAR + + # If no shock conditions, delegate to the reduced-form method + if shock_conditions is None: + fitted = FittedVAR.model_construct( + idata=self.idata, + n_lags=self.n_lags, + data=self.data, + var_names=self.var_names, + ) + return fitted.conditional_forecast(steps=steps, conditions=conditions, exog_future=exog_future) + + if exog_future is not None: + raise NotImplementedError( + "exog_future is not yet supported with shock_conditions. " + "Use shock_conditions=None to use the reduced-form path with exog support." + ) + + return self._structural_conditional_forecast(steps, conditions, shock_conditions) + + def _validate_structural_conditions( + self, conditions: "list[ForecastCondition]", shock_conditions: "list[ForecastCondition]", steps: int + ) -> list[str]: + """Validate observable and shock conditions, returning shock names. + + Args: + conditions: Observable variable conditions. + shock_conditions: Structural shock conditions. + steps: Number of forecast steps. + + Returns: + List of shock variable names from the structural matrix. + """ + for cond in conditions: + if cond.variable not in self.var_names: + raise ValueError(f"Condition variable '{cond.variable}' not in var_names") + for p in cond.periods: + if p < 0 or p >= steps: + raise ValueError(f"Condition period {p} out of range for {steps} steps") + + shock_names = self.idata.posterior["structural_shock_matrix"].coords["shock"].values.tolist() + for cond in shock_conditions: + if cond.variable not in shock_names: + raise ValueError(f"Shock condition variable '{cond.variable}' not in shock_names {shock_names}") + return shock_names + + @staticmethod + def _build_structural_constraint_system( + conditions: "list[ForecastCondition]", + shock_conditions: "list[ForecastCondition]", + var_names: list[str], + shock_names: list[str], + ma_coefficients: list[np.ndarray], + P: np.ndarray, + unconditional: np.ndarray, + steps: int, + n_vars: int, + ) -> tuple[np.ndarray, np.ndarray]: + """Build the combined observable + shock constraint system. + + Returns: + Tuple of (R, target) arrays for the constraint system. + """ + constraint_rows = [] + constraint_targets = [] + + # Observable constraints + for cond in conditions: + var_idx = var_names.index(cond.variable) + for period, value in zip(cond.periods, cond.values, strict=True): + row = np.zeros(steps * n_vars) + for s in range(period + 1): + structural_response = ma_coefficients[period - s] @ P + row[s * n_vars : (s + 1) * n_vars] = structural_response[var_idx, :] + constraint_rows.append(row) + constraint_targets.append(value - unconditional[period, var_idx]) + + # Shock constraints + for cond in shock_conditions: + shock_idx = shock_names.index(cond.variable) + for period, value in zip(cond.periods, cond.values, strict=True): + row = np.zeros(steps * n_vars) + row[period * n_vars + shock_idx] = 1.0 + constraint_rows.append(row) + constraint_targets.append(value) + + return np.array(constraint_rows), np.array(constraint_targets) + + def _structural_conditional_forecast( + self, steps: int, conditions: "list[ForecastCondition]", shock_conditions: "list[ForecastCondition]" + ) -> "ConditionalForecastResult": + """Compute structural conditional forecast with shock constraints.""" + from impulso.fitted import FittedVAR + from impulso.results import ConditionalForecastResult + + shock_names = self._validate_structural_conditions(conditions, shock_conditions, steps) + + B_draws = self.idata.posterior["B"].values + intercept_draws = self.idata.posterior["intercept"].values + P_draws = self.idata.posterior["structural_shock_matrix"].values + n_chains, n_draws, n_vars, _ = B_draws.shape + + y_hist = self.data.endog[-self.n_lags :] + forecasts = np.zeros((n_chains, n_draws, steps, n_vars)) + + for c in range(n_chains): + for d in range(n_draws): + B = B_draws[c, d] + P = P_draws[c, d] + + ma_coefficients = FittedVAR._ma_coefficients_single(B, n_vars, self.n_lags, steps) + unconditional = self._unconditional_forecast_single( + B, intercept_draws[c, d], y_hist, steps, self.n_lags, n_vars + ) + + R, target = self._build_structural_constraint_system( + conditions, + shock_conditions, + self.var_names, + shock_names, + ma_coefficients, + P, + unconditional, + steps, + n_vars, + ) + + structural_shocks, _, _, _ = np.linalg.lstsq(R, target, rcond=None) + structural_shocks = structural_shocks.reshape(steps, n_vars) + + conditional = unconditional.copy() + for h in range(steps): + for s in range(h + 1): + conditional[h] += ma_coefficients[h - s] @ P @ structural_shocks[s] + + forecasts[c, d] = conditional + + forecast_da = xr.DataArray( + forecasts, + dims=["chain", "draw", "step", "variable"], + coords={"variable": self.var_names}, + name="forecast", + ) + idata = az.InferenceData(posterior_predictive=xr.Dataset({"forecast": forecast_da})) + all_conditions = conditions + (shock_conditions or []) + return ConditionalForecastResult(idata=idata, steps=steps, var_names=self.var_names, conditions=all_conditions) + + @staticmethod + def _unconditional_forecast_single( + B: np.ndarray, + intercept: np.ndarray, + y_hist: np.ndarray, + steps: int, + n_lags: int, + n_vars: int, + ) -> np.ndarray: + """Compute unconditional forecast for a single posterior draw. + + Returns: + Array of shape (steps, n_vars). + """ + y_buffer = y_hist.copy() + unconditional = np.zeros((steps, n_vars)) + for h in range(steps): + x_lag = np.concatenate([y_buffer[-(lag + 1)] for lag in range(n_lags)]) + unconditional[h] = intercept + B @ x_lag + y_buffer = np.vstack([y_buffer[1:], unconditional[h].reshape(1, -1)]) + return unconditional diff --git a/src/impulso/results.py b/src/impulso/results.py index ba2a9da..97f74ff 100644 --- a/src/impulso/results.py +++ b/src/impulso/results.py @@ -101,6 +101,16 @@ def plot(self) -> Figure: return plot_forecast(self) +class ConditionalForecastResult(ForecastResult): + """Result from conditional VAR forecasting. + + Attributes: + conditions: List of ForecastConditions applied. + """ + + conditions: list # list[ForecastCondition] — bare list avoids Pydantic rebuild issues + + class IRFResult(VARResultBase): """Result from impulse response function computation. diff --git a/tests/test_conditional_forecast.py b/tests/test_conditional_forecast.py new file mode 100644 index 0000000..7329846 --- /dev/null +++ b/tests/test_conditional_forecast.py @@ -0,0 +1,126 @@ +"""Tests for conditional forecasting.""" + +import numpy as np +import pytest +from pydantic import ValidationError + +from impulso.conditions import ForecastCondition + + +class TestForecastCondition: + def test_basic_construction(self): + fc = ForecastCondition(variable="y1", periods=[0, 1, 2], values=[1.0, 1.0, 1.0]) + assert fc.variable == "y1" + assert fc.periods == [0, 1, 2] + assert fc.constraint_type == "hard" + + def test_frozen(self): + fc = ForecastCondition(variable="y1", periods=[0], values=[1.0]) + with pytest.raises(ValidationError): + fc.variable = "y2" + + def test_rejects_mismatched_lengths(self): + with pytest.raises(ValidationError, match="periods length"): + ForecastCondition(variable="y1", periods=[0, 1], values=[1.0]) + + def test_rejects_empty_periods(self): + with pytest.raises(ValidationError, match="non-empty"): + ForecastCondition(variable="y1", periods=[], values=[]) + + def test_rejects_negative_periods(self): + with pytest.raises(ValidationError, match="non-negative"): + ForecastCondition(variable="y1", periods=[-1], values=[1.0]) + + +class TestConditionalForecastOnFittedVAR: + def test_returns_forecast_result(self, var_data_2v): + from impulso.conjugate import ConjugateVAR + + cvar = ConjugateVAR(lags=1, draws=50, random_seed=42) + fitted = cvar.fit(var_data_2v) + conditions = [ + ForecastCondition(variable="y1", periods=[0, 1, 2, 3], values=[0.5, 0.5, 0.5, 0.5]), + ] + result = fitted.conditional_forecast(steps=8, conditions=conditions) + assert result.median().shape == (8, 2) + + def test_constrained_periods_match_target(self, var_data_2v): + from impulso.conjugate import ConjugateVAR + + target = 0.5 + cvar = ConjugateVAR(lags=1, draws=50, random_seed=42) + fitted = cvar.fit(var_data_2v) + conditions = [ + ForecastCondition(variable="y1", periods=[0, 1], values=[target, target]), + ] + result = fitted.conditional_forecast(steps=4, conditions=conditions) + median = result.median() + # Constrained periods should be close to target (exact for hard constraints) + np.testing.assert_allclose(median.iloc[0]["y1"], target, atol=1e-6) + np.testing.assert_allclose(median.iloc[1]["y1"], target, atol=1e-6) + + def test_unconstrained_variable_differs_from_unconditional(self, var_data_2v): + from impulso.conjugate import ConjugateVAR + + cvar = ConjugateVAR(lags=1, draws=50, random_seed=42) + fitted = cvar.fit(var_data_2v) + conditions = [ + ForecastCondition(variable="y1", periods=[0, 1, 2, 3], values=[5.0, 5.0, 5.0, 5.0]), + ] + unconditional = fitted.forecast(steps=4).median() + conditional = fitted.conditional_forecast(steps=4, conditions=conditions).median() + # y2 should differ because y1 is forced far from its unconditional path + assert not np.allclose(unconditional["y2"].values, conditional["y2"].values) + + def test_rejects_unknown_variable(self, var_data_2v): + from impulso.conjugate import ConjugateVAR + + cvar = ConjugateVAR(lags=1, draws=50, random_seed=42) + fitted = cvar.fit(var_data_2v) + conditions = [ + ForecastCondition(variable="unknown", periods=[0], values=[1.0]), + ] + with pytest.raises(ValueError, match="unknown"): + fitted.conditional_forecast(steps=4, conditions=conditions) + + def test_rejects_period_out_of_range(self, var_data_2v): + from impulso.conjugate import ConjugateVAR + + cvar = ConjugateVAR(lags=1, draws=50, random_seed=42) + fitted = cvar.fit(var_data_2v) + conditions = [ + ForecastCondition(variable="y1", periods=[10], values=[1.0]), + ] + with pytest.raises(ValueError, match="out of range"): + fitted.conditional_forecast(steps=4, conditions=conditions) + + +class TestConditionalForecastOnIdentifiedVAR: + def test_returns_forecast_result(self, var_data_2v): + from impulso.conjugate import ConjugateVAR + from impulso.identification import Cholesky + + cvar = ConjugateVAR(lags=1, draws=50, random_seed=42) + fitted = cvar.fit(var_data_2v) + identified = fitted.set_identification_strategy(Cholesky(ordering=["y1", "y2"])) + conditions = [ + ForecastCondition(variable="y1", periods=[0, 1], values=[0.5, 0.5]), + ] + result = identified.conditional_forecast(steps=4, conditions=conditions) + assert result.median().shape == (4, 2) + + def test_with_shock_conditions(self, var_data_2v): + from impulso.conjugate import ConjugateVAR + from impulso.identification import Cholesky + + cvar = ConjugateVAR(lags=1, draws=50, random_seed=42) + fitted = cvar.fit(var_data_2v) + identified = fitted.set_identification_strategy(Cholesky(ordering=["y1", "y2"])) + conditions = [ + ForecastCondition(variable="y1", periods=[0, 1], values=[0.5, 0.5]), + ] + shock_conditions = [ + ForecastCondition(variable="y1", periods=[0, 1], values=[0.0, 0.0]), + ] + result = identified.conditional_forecast(steps=4, conditions=conditions, shock_conditions=shock_conditions) + assert result.median().shape == (4, 2) diff --git a/tests/test_conjugate.py b/tests/test_conjugate.py new file mode 100644 index 0000000..d324287 --- /dev/null +++ b/tests/test_conjugate.py @@ -0,0 +1,145 @@ +"""Tests for ConjugateVAR (direct NIW posterior sampling).""" + +import numpy as np +import pytest +from pydantic import ValidationError + +from impulso.priors import MinnesotaPrior + + +class TestConjugateVARConstruction: + def test_basic_construction(self): + from impulso.conjugate import ConjugateVAR + + cvar = ConjugateVAR(lags=2) + assert cvar.lags == 2 + assert cvar.draws == 2000 + + def test_custom_prior(self): + from impulso.conjugate import ConjugateVAR + + prior = MinnesotaPrior(tightness=0.2, cross_shrinkage=0.3) + cvar = ConjugateVAR(lags=2, prior=prior) + assert cvar.prior == prior + + def test_frozen(self): + from impulso.conjugate import ConjugateVAR + + cvar = ConjugateVAR(lags=2) + with pytest.raises(ValidationError): + cvar.lags = 4 + + def test_rejects_negative_draws(self): + from impulso.conjugate import ConjugateVAR + + with pytest.raises(ValidationError): + ConjugateVAR(lags=2, draws=0) + + +class TestConjugateVARFit: + def test_fit_returns_fitted_var(self, var_data_2v): + from impulso.conjugate import ConjugateVAR + from impulso.fitted import FittedVAR + + cvar = ConjugateVAR(lags=1, draws=100, random_seed=42) + fitted = cvar.fit(var_data_2v) + assert isinstance(fitted, FittedVAR) + + def test_idata_has_required_variables(self, var_data_2v): + from impulso.conjugate import ConjugateVAR + + cvar = ConjugateVAR(lags=1, draws=100, random_seed=42) + fitted = cvar.fit(var_data_2v) + assert "B" in fitted.idata.posterior + assert "intercept" in fitted.idata.posterior + assert "Sigma" in fitted.idata.posterior + + def test_posterior_shapes(self, var_data_2v): + from impulso.conjugate import ConjugateVAR + + n_draws = 100 + cvar = ConjugateVAR(lags=2, draws=n_draws, random_seed=42) + fitted = cvar.fit(var_data_2v) + B = fitted.idata.posterior["B"].values + assert B.shape == (1, n_draws, 2, 4) + intercept = fitted.idata.posterior["intercept"].values + assert intercept.shape == (1, n_draws, 2) + sigma = fitted.idata.posterior["Sigma"].values + assert sigma.shape == (1, n_draws, 2, 2) + + def test_sigma_positive_definite(self, var_data_2v): + from impulso.conjugate import ConjugateVAR + + cvar = ConjugateVAR(lags=1, draws=100, random_seed=42) + fitted = cvar.fit(var_data_2v) + sigma = fitted.idata.posterior["Sigma"].values + for d in range(sigma.shape[1]): + eigvals = np.linalg.eigvalsh(sigma[0, d]) + assert np.all(eigvals > 0), f"Draw {d} has non-positive eigenvalue" + + def test_sigma_symmetric(self, var_data_2v): + from impulso.conjugate import ConjugateVAR + + cvar = ConjugateVAR(lags=1, draws=100, random_seed=42) + fitted = cvar.fit(var_data_2v) + sigma = fitted.idata.posterior["Sigma"].values + np.testing.assert_allclose(sigma, np.swapaxes(sigma, -2, -1), atol=1e-10) + + def test_reproducible_with_seed(self, var_data_2v): + from impulso.conjugate import ConjugateVAR + + cvar1 = ConjugateVAR(lags=1, draws=50, random_seed=123) + cvar2 = ConjugateVAR(lags=1, draws=50, random_seed=123) + fitted1 = cvar1.fit(var_data_2v) + fitted2 = cvar2.fit(var_data_2v) + np.testing.assert_array_equal( + fitted1.idata.posterior["B"].values, + fitted2.idata.posterior["B"].values, + ) + + def test_var_names_correct(self, var_data_2v): + from impulso.conjugate import ConjugateVAR + + cvar = ConjugateVAR(lags=1, draws=50, random_seed=42) + fitted = cvar.fit(var_data_2v) + assert fitted.var_names == ["y1", "y2"] + + def test_n_lags_stored(self, var_data_2v): + from impulso.conjugate import ConjugateVAR + + cvar = ConjugateVAR(lags=3, draws=50, random_seed=42) + fitted = cvar.fit(var_data_2v) + assert fitted.n_lags == 3 + + def test_downstream_forecast_works(self, var_data_2v): + from impulso.conjugate import ConjugateVAR + + cvar = ConjugateVAR(lags=1, draws=50, random_seed=42) + fitted = cvar.fit(var_data_2v) + result = fitted.forecast(steps=4) + assert result.median().shape == (4, 2) + + def test_downstream_identification_works(self, var_data_2v): + from impulso.conjugate import ConjugateVAR + from impulso.identification import Cholesky + + cvar = ConjugateVAR(lags=1, draws=50, random_seed=42) + fitted = cvar.fit(var_data_2v) + identified = fitted.set_identification_strategy(Cholesky(ordering=["y1", "y2"])) + irfs = identified.impulse_response(horizon=10) + assert irfs.median().shape[0] == 11 + + def test_lag_selection_string(self, var_data_2v): + from impulso.conjugate import ConjugateVAR + + cvar = ConjugateVAR(lags="bic", draws=50, random_seed=42) + fitted = cvar.fit(var_data_2v) + assert fitted.n_lags >= 1 + + def test_works_with_dummy_observations(self, var_data_2v): + from impulso.conjugate import ConjugateVAR + + augmented = var_data_2v.with_dummy_observations(n_lags=2, mu=5.0, delta=1.0) + cvar = ConjugateVAR(lags=2, draws=50, random_seed=42) + fitted = cvar.fit(augmented) + assert fitted.idata.posterior["B"].values.shape == (1, 50, 2, 4) diff --git a/tests/test_dummy_observations.py b/tests/test_dummy_observations.py new file mode 100644 index 0000000..a62aa81 --- /dev/null +++ b/tests/test_dummy_observations.py @@ -0,0 +1,93 @@ +"""Tests for dummy observation priors on VARData.""" + +import numpy as np +import pandas as pd +import pytest + +from impulso.data import VARData + + +@pytest.fixture +def var_data(): + rng = np.random.default_rng(42) + T, n = 100, 3 + endog = rng.standard_normal((T, n)) + index = pd.date_range("2000-01-01", periods=T, freq="QS") + return VARData(endog=endog, endog_names=["gdp", "inflation", "rate"], index=index) + + +class TestDummyObservationPriors: + def test_sum_of_coefficients_appends_n_vars_rows(self, var_data): + augmented = var_data.with_dummy_observations(n_lags=4, mu=5.0) + assert augmented.endog.shape[0] == var_data.endog.shape[0] + 3 + + def test_single_unit_root_appends_one_row(self, var_data): + augmented = var_data.with_dummy_observations(n_lags=4, delta=1.0) + assert augmented.endog.shape[0] == var_data.endog.shape[0] + 1 + + def test_both_dummies_append_correct_rows(self, var_data): + augmented = var_data.with_dummy_observations(n_lags=4, mu=5.0, delta=1.0) + assert augmented.endog.shape[0] == var_data.endog.shape[0] + 4 + + def test_sum_of_coefficients_values(self, var_data): + mu = 5.0 + augmented = var_data.with_dummy_observations(n_lags=4, mu=mu) + y_bar = var_data.endog.mean(axis=0) + dummy_rows = augmented.endog[var_data.endog.shape[0] :] + for i in range(3): + expected = np.zeros(3) + expected[i] = y_bar[i] / mu + np.testing.assert_allclose(dummy_rows[i], expected) + + def test_single_unit_root_values(self, var_data): + delta = 1.0 + augmented = var_data.with_dummy_observations(n_lags=4, delta=delta) + y_bar = var_data.endog.mean(axis=0) + dummy_row = augmented.endog[-1] + np.testing.assert_allclose(dummy_row, y_bar / delta) + + def test_preserves_original_data(self, var_data): + augmented = var_data.with_dummy_observations(n_lags=4, mu=5.0) + np.testing.assert_array_equal(augmented.endog[: var_data.endog.shape[0]], var_data.endog) + + def test_returns_new_vardata(self, var_data): + augmented = var_data.with_dummy_observations(n_lags=4, mu=5.0) + assert augmented is not var_data + assert isinstance(augmented, VARData) + + def test_index_extended(self, var_data): + augmented = var_data.with_dummy_observations(n_lags=4, mu=5.0) + assert len(augmented.index) == augmented.endog.shape[0] + + def test_endog_names_preserved(self, var_data): + augmented = var_data.with_dummy_observations(n_lags=4, mu=5.0) + assert augmented.endog_names == var_data.endog_names + + def test_raises_if_neither_mu_nor_delta(self, var_data): + with pytest.raises(ValueError, match="At least one"): + var_data.with_dummy_observations(n_lags=4) + + def test_raises_if_mu_not_positive(self, var_data): + with pytest.raises(ValueError, match="mu must be"): + var_data.with_dummy_observations(n_lags=4, mu=-1.0) + + def test_raises_if_delta_not_positive(self, var_data): + with pytest.raises(ValueError, match="delta must be"): + var_data.with_dummy_observations(n_lags=4, delta=0.0) + + def test_raises_if_n_lags_not_positive(self, var_data): + with pytest.raises(ValueError, match="n_lags must be"): + var_data.with_dummy_observations(n_lags=0, mu=5.0) + + def test_exog_zero_padded(self): + rng = np.random.default_rng(42) + T, n = 100, 2 + endog = rng.standard_normal((T, n)) + exog = rng.standard_normal((T, 1)) + index = pd.date_range("2000-01-01", periods=T, freq="QS") + data = VARData(endog=endog, endog_names=["y1", "y2"], exog=exog, exog_names=["x1"], index=index) + augmented = data.with_dummy_observations(n_lags=2, mu=5.0) + assert augmented.exog.shape[0] == augmented.endog.shape[0] + np.testing.assert_array_equal(augmented.exog[:T], data.exog) + np.testing.assert_array_equal(augmented.exog[T:], 0.0) + assert augmented.exog_names == ["x1"] diff --git a/tests/test_glp.py b/tests/test_glp.py new file mode 100644 index 0000000..16e4e85 --- /dev/null +++ b/tests/test_glp.py @@ -0,0 +1,82 @@ +"""Tests for GLP hierarchical prior selection on ConjugateVAR.""" + +import numpy as np + +from impulso.priors import MinnesotaPrior + + +class TestMarginalLikelihood: + def test_returns_finite_scalar(self, var_data_2v): + from impulso.conjugate import ConjugateVAR + + cvar = ConjugateVAR(lags=1) + ml = cvar.marginal_likelihood(var_data_2v) + assert np.isfinite(ml) + + def test_varies_with_prior(self, var_data_2v): + from impulso.conjugate import ConjugateVAR + + cvar_tight = ConjugateVAR(lags=1, prior=MinnesotaPrior(tightness=0.01)) + cvar_loose = ConjugateVAR(lags=1, prior=MinnesotaPrior(tightness=1.0)) + ml_tight = cvar_tight.marginal_likelihood(var_data_2v) + ml_loose = cvar_loose.marginal_likelihood(var_data_2v) + assert ml_tight != ml_loose + + def test_higher_for_true_lag_order(self, var_data_2v): + """Marginal likelihood should favour the true DGP lag order (1).""" + from impulso.conjugate import ConjugateVAR + + ml_1 = ConjugateVAR(lags=1).marginal_likelihood(var_data_2v) + ml_8 = ConjugateVAR(lags=8).marginal_likelihood(var_data_2v) + assert ml_1 > ml_8 + + +class TestOptimizePrior: + def test_returns_minnesota_prior(self, var_data_2v): + from impulso.conjugate import ConjugateVAR + + cvar = ConjugateVAR(lags=1) + optimal = cvar.optimize_prior(var_data_2v) + assert isinstance(optimal, MinnesotaPrior) + + def test_optimal_tightness_positive(self, var_data_2v): + from impulso.conjugate import ConjugateVAR + + cvar = ConjugateVAR(lags=1) + optimal = cvar.optimize_prior(var_data_2v) + assert optimal.tightness > 0 + + def test_optimal_cross_shrinkage_in_bounds(self, var_data_2v): + from impulso.conjugate import ConjugateVAR + + cvar = ConjugateVAR(lags=1) + optimal = cvar.optimize_prior(var_data_2v) + assert 0.01 <= optimal.cross_shrinkage <= 1.0 + + def test_optimal_has_higher_marginal_likelihood(self, var_data_2v): + from impulso.conjugate import ConjugateVAR + + cvar_default = ConjugateVAR(lags=1) + ml_default = cvar_default.marginal_likelihood(var_data_2v) + + optimal_prior = cvar_default.optimize_prior(var_data_2v) + cvar_optimal = ConjugateVAR(lags=1, prior=optimal_prior) + ml_optimal = cvar_optimal.marginal_likelihood(var_data_2v) + + assert ml_optimal >= ml_default - 1e-6 # allow tiny numerical tolerance + + def test_preserves_decay_setting(self, var_data_2v): + from impulso.conjugate import ConjugateVAR + + cvar = ConjugateVAR(lags=1, prior=MinnesotaPrior(decay="geometric")) + optimal = cvar.optimize_prior(var_data_2v) + assert optimal.decay == "geometric" + + def test_minnesota_optimized_shorthand(self, var_data_2v): + """prior='minnesota_optimized' should trigger automatic optimisation.""" + from impulso.conjugate import ConjugateVAR + from impulso.fitted import FittedVAR + + cvar = ConjugateVAR(lags=1, prior="minnesota_optimized", draws=50, random_seed=42) + fitted = cvar.fit(var_data_2v) + assert isinstance(fitted, FittedVAR) diff --git a/tests/test_long_run_restriction.py b/tests/test_long_run_restriction.py new file mode 100644 index 0000000..d2f3653 --- /dev/null +++ b/tests/test_long_run_restriction.py @@ -0,0 +1,150 @@ +"""Tests for Blanchard-Quah long-run identification.""" + +import arviz as az +import numpy as np +import pytest +import xarray as xr +from pydantic import ValidationError + +from impulso.protocols import IdentificationScheme + + +@pytest.fixture +def stationary_idata_2v(): + """Synthetic InferenceData with stationary VAR(1) coefficients.""" + rng = np.random.default_rng(42) + n_chains, n_draws, n_vars = 2, 50, 2 + + # Stationary coefficients: eigenvalues inside unit circle + B = np.zeros((n_chains, n_draws, n_vars, n_vars)) + for c in range(n_chains): + for d in range(n_draws): + # Diagonal with small values ensures stationarity + B[c, d] = np.diag(rng.uniform(0.1, 0.4, n_vars)) + + intercept = rng.standard_normal((n_chains, n_draws, n_vars)) * 0.01 + sigma = np.zeros((n_chains, n_draws, n_vars, n_vars)) + for c in range(n_chains): + for d in range(n_draws): + A = rng.standard_normal((n_vars, n_vars)) * 0.5 + sigma[c, d] = A @ A.T + np.eye(n_vars) + + posterior = xr.Dataset({ + "B": xr.DataArray(B, dims=["chain", "draw", "var", "coeff"]), + "intercept": xr.DataArray(intercept, dims=["chain", "draw", "var"]), + "Sigma": xr.DataArray(sigma, dims=["chain", "draw", "var1", "var2"]), + }) + return az.InferenceData(posterior=posterior) + + +class TestLongRunRestrictionConstruction: + def test_basic_construction(self): + from impulso.identification import LongRunRestriction + + lr = LongRunRestriction(ordering=["output", "prices"]) + assert lr.ordering == ["output", "prices"] + + def test_frozen(self): + from impulso.identification import LongRunRestriction + + lr = LongRunRestriction(ordering=["a", "b"]) + with pytest.raises(ValidationError): + lr.ordering = ["b", "a"] + + def test_satisfies_protocol(self): + from impulso.identification import LongRunRestriction + + lr = LongRunRestriction(ordering=["a", "b"]) + assert isinstance(lr, IdentificationScheme) + + +class TestLongRunRestrictionIdentify: + def test_produces_structural_shock_matrix(self, stationary_idata_2v): + from impulso.identification import LongRunRestriction + + lr = LongRunRestriction(ordering=["y1", "y2"]) + result = lr.identify(stationary_idata_2v, var_names=["y1", "y2"]) + assert "structural_shock_matrix" in result.posterior + + def test_output_shape(self, stationary_idata_2v): + from impulso.identification import LongRunRestriction + + lr = LongRunRestriction(ordering=["y1", "y2"]) + result = lr.identify(stationary_idata_2v, var_names=["y1", "y2"]) + P = result.posterior["structural_shock_matrix"].values + assert P.shape == (2, 50, 2, 2) + + def test_no_nan_values(self, stationary_idata_2v): + from impulso.identification import LongRunRestriction + + lr = LongRunRestriction(ordering=["y1", "y2"]) + result = lr.identify(stationary_idata_2v, var_names=["y1", "y2"]) + P = result.posterior["structural_shock_matrix"].values + assert not np.any(np.isnan(P)) + + def test_long_run_impact_is_lower_triangular(self, stationary_idata_2v): + """The long-run cumulative impact C(1) @ P should be lower triangular.""" + from impulso.identification import LongRunRestriction + + lr = LongRunRestriction(ordering=["y1", "y2"]) + result = lr.identify(stationary_idata_2v, var_names=["y1", "y2"]) + P = result.posterior["structural_shock_matrix"].values + B = stationary_idata_2v.posterior["B"].values + n_vars = 2 + + for c in range(2): + for d in range(50): + lag_coefficient_sum = B[c, d, :, :n_vars] + long_run_multiplier = np.linalg.inv(np.eye(n_vars) - lag_coefficient_sum) + long_run_impact = long_run_multiplier @ P[c, d] + # Upper triangle (excluding diagonal) should be ~zero + np.testing.assert_allclose( + np.triu(long_run_impact, k=1), + 0.0, + atol=1e-10, + ) + + def test_reordering_works(self, stationary_idata_2v): + from impulso.identification import LongRunRestriction + + lr = LongRunRestriction(ordering=["y2", "y1"]) + result = lr.identify(stationary_idata_2v, var_names=["y1", "y2"]) + assert result.posterior["structural_shock_matrix"].coords["shock"].values.tolist() == ["y2", "y1"] + + def test_coordinates_match_ordering(self, stationary_idata_2v): + from impulso.identification import LongRunRestriction + + lr = LongRunRestriction(ordering=["y1", "y2"]) + result = lr.identify(stationary_idata_2v, var_names=["y1", "y2"]) + assert result.posterior["structural_shock_matrix"].coords["shock"].values.tolist() == ["y1", "y2"] + assert result.posterior["structural_shock_matrix"].coords["response"].values.tolist() == ["y1", "y2"] + + def test_preserves_other_posterior_variables(self, stationary_idata_2v): + from impulso.identification import LongRunRestriction + + lr = LongRunRestriction(ordering=["y1", "y2"]) + result = lr.identify(stationary_idata_2v, var_names=["y1", "y2"]) + assert "B" in result.posterior + assert "Sigma" in result.posterior + + def test_raises_when_B_missing(self, stationary_idata_2v): + from impulso.identification import LongRunRestriction + + lr = LongRunRestriction(ordering=["y1", "y2"]) + idata_no_B = az.InferenceData(posterior=stationary_idata_2v.posterior.drop_vars("B")) + with pytest.raises(ValueError, match="requires 'B'"): + lr.identify(idata_no_B, var_names=["y1", "y2"]) + + def test_raises_for_unknown_variables(self, stationary_idata_2v): + from impulso.identification import LongRunRestriction + + lr = LongRunRestriction(ordering=["y1", "unknown"]) + with pytest.raises(ValueError, match="unknown variables"): + lr.identify(stationary_idata_2v, var_names=["y1", "y2"]) + + def test_raises_for_wrong_ordering_length(self, stationary_idata_2v): + from impulso.identification import LongRunRestriction + + lr = LongRunRestriction(ordering=["y1"]) + with pytest.raises(ValueError, match="exactly 2 variables"): + lr.identify(stationary_idata_2v, var_names=["y1", "y2"])