Skewed Constraints for Conditional Forecasting¶

Paul Labonne, Andrea Renzetti

Introduction¶

This document describes conditional forecasting with skewed (non-normal) constraint distributions in the BVAR model. The approach lets users impose distributional constraints on selected variables while accounting for parameter and shock uncertainty.

Setup¶

The conditional forecast is constructed as:

where: - \(y\) is the \((Hn \times 1)\) forecast vector over \(H\) steps - \(b\) is the \((Hn \times 1)\) deterministic component (from constant and lagged values) - \(M\) is the \((Hn \times Hn)\) block lower-triangular impulse response matrix - \(\varepsilon\) is the \((Hn \times 1)\) vector of structural shocks

We impose distributional constraints:

where: - \(C\) is the \((m \times Hn)\) selection matrix (m = number of constraints) - \(SN_{\text{DP}}(\cdot)\) is the Skew Normal distribution (Azzalini & Capitanio, 1999, Section 4.1) - \(f\) is the location vector (constraint mean) - \(\Sigma_f\) is the scale matrix (constraint covariance) - \(\alpha\) is the shape vector (constraint skewness)

Algorithm: Shock Adjustment with Constraint Sampling¶

This is the algorithm forecast() runs by default (method="andersson_et_al"), following Algorithm 1 of Andersson, Palmqvist & Waggoner (2010).

Choosing method¶

Key Idea¶

We decompose the structural shocks into an unconstrained component and a constraint-satisfying adjustment:

where \(\xi \sim N(0, I_{Hn})\) represents unconstrained shocks and \(\psi\) is the adjustment ensuring \(Cy = z\) exactly.

Four-Step Algorithm¶

Step 1: Draw unconstrained shocks

Step 2: Draw constraint targets from the constraint distribution

where \(z\) represents the target value for \(Cy - Cb = D\varepsilon\) (with \(D = CM^T\)).

Step 3: Adjust shocks to satisfy constraints

where \(D^{*} = D^T (D D^T)^{-1}\) is the Moore-Penrose pseudo-inverse of the wide, full-row-rank matrix \(D\) (computed via SVD, numpy.linalg.pinv).

Step 4: Compute the forecast

Mathematical Intuition¶

Starting from the constraint requirement \(Cy = z\):

We decompose \(\varepsilon = \xi + \psi\) and solve for \(\psi\):

Thus:

Key property: \((I - D^{*}D)\) projects onto the null space of \(D\), so unconstrained directions in \(\xi\) only affect the forecast through shock directions orthogonal to the constraints.

Implementation Paths¶

Path 1: Conditional Mean Only (point_only=True)¶

Computes the conditional forecast mean without sampling:

Use this for fast point estimates when uncertainty quantification is not needed.

Path 2: Default Constraint Sampling¶

If no custom sampler is provided, constraint targets are drawn from a Skew Normal distribution:

When all shape parameters satisfy \(\alpha_i = 0\), the distribution becomes Normal. The method supports independent distributions for each constraint.

Path 3: Custom Constraint Sampler¶

Users can provide a constraint_sampler callable that returns draws from an arbitrary constraint distribution:

def my_constraint_sampler():
    # User-defined sampling logic
    return constraint_draw  # Shape: (n_constraints,)


forecast = bvar.forecast(
    H=8,
    constraint_mean=constraint_mean,
    constraint_variance=constraint_variance,
    constraint_sampler=my_constraint_sampler,
)

This enables flexible specifications such as truncated distributions or mixtures, with the constraint values automatically transformed internally to \(z - Cb\) for the shock constraint.

Constraint Specification¶

Constraints are provided to the forecast() method as three \((H \times n)\) arrays:

constraint_mean = np.full((H, n), np.nan)  # NaNs for unconstrained
constraint_mean[0:2, 0] = 2.5  # Fix variable 0 to 2.5 for steps 0–1

constraint_variance = np.full((H, n), np.nan)
constraint_variance[:, 0] = 0.5  # Allow ±0.71 standard deviation for variable 0

constraint_shape = np.full((H, n), np.nan)
constraint_shape[:, 0] = 2.0  # Right-skewed constraint for variable 0

The get_constraint() function internally converts these arrays into: - Selection matrix \(C\) (shape: \(m \times Hn\)) - Location vector \(f\) (shape: \(m\)) - Scale matrix \(\Sigma_f\) (shape: \(m \times m\), diagonal) - Shape vector \(\alpha\) (shape: \(m\))

where \(m\) is the number of non-NaN constraint entries.

Gibbs Sampler for Parameter Uncertainty¶

To account for uncertainty in both the VAR parameters and the constraint distribution, the conditional forecasting procedure embeds the shock adjustment algorithm within a Gibbs sampler:

Initialization: Use posterior point estimates \(\beta_{\text{point}}\) and \(\Sigma_{\text{point}}\).

For each MCMC iteration \(i = 1, \ldots, N_{\text{draws}}\):

1. Draw constrained forecast \(y^{(i)}\) using the shock adjustment algorithm above
2. Augment the dataset with the conditional forecast: \(y^{(i)}_{\text{data}} = [y_{\text{obs}}, y^{(i)}]\)
3. Re-estimate the BVAR on the augmented data using its posterior sampler with the full state to obtain posterior draws \((\beta^{(i)}, \Sigma^{(i)})\)
4. Use \((\beta^{(i)}, \Sigma^{(i)})\) for the next iteration

Post-processing: Discard the first \(N_{\text{burn}}\) draws (default: \(N_{\text{draws}} / 2\)).

This ensures internal consistency: the conditional forecast paths follow the BVAR dynamics while propagating parameter and forecast uncertainty.

Alternative Algorithms (Reference)¶

Method 1: Affine Transformation of Skew Normal (Antolín-Díaz et al., 2021)¶

If \(X \sim SN_{\text{DP}}(\mu, \Omega, \alpha)\), then for any matrix \(A\) and vector \(c\):

where the shape parameter transforms according to equation (10) of Azzalini & Capitanio (1999).

Drawback: Computing the posterior shock covariance \(\Omega_\varepsilon = D^{*}\Sigma_f(D^{*})^T + I - D^{*}DD^T(D^{*})^T\) can produce non-positive-definite matrices when variance constraints are very tight, causing numerical instability.

Reference algorithm: explicit null-space decomposition¶

Decomposes the shock vector into constrained and unconstrained components using the null space of \(D\):

where \(\hat{D}\) is an orthonormal basis for \(\text{null}(D)\).

Drawback: Requires explicit null space computation at each iteration, which is slower than the shock adjustment method.

Usage Example¶

import numpy as np
import bvar as bv

# Estimate model
bvar = bv.BVAR(
    n_lags=2,
    model=bv.NaturalConjugate(minnesota=True, soc=True, sur=True),
    stationary=False,
    random_state=42,
)
# default hyperparameters (no optimise_hyperparameters call) to keep the example short
bvar.sample(data, N_draws=5000)

# Set up constraints: fix GDP growth at 2% for next 2 quarters
H = 8
constraint_mean = np.full((H, bvar.n), np.nan)
constraint_mean[0:2, 0] = 2.0  # GDP index

constraint_variance = np.full((H, bvar.n), np.nan)
constraint_variance[0:2, 0] = 0.25  # Tight around the mean

constraint_shape = np.full((H, bvar.n), np.nan)
constraint_shape[0:2, 0] = 0.0  # Normal distribution (no skew)

# Generate conditional forecast
bvar.forecast(
    H=H,
    constraint_mean=constraint_mean,
    constraint_variance=constraint_variance,
    constraint_shape=constraint_shape,
    method="andersson_et_al",
    N_draws=5000,
    progressbar=True,
)

# Extract results
y_cond = (
    bvar.forecast_conditional
)  # Shape: (min(N_draws, self.N_draws) - N_burn, T+H, n)

References¶

• Andersson, M. K., Palmqvist, S., & Waggoner, D. F. (2010). Density-conditional forecasts in dynamic multivariate models. Sveriges Riksbank Working Paper Series, 243.
• Antolín-Díaz, J., Petrella, I., & Rubio-Ramírez, J. F. (2021). Structural scenario analysis with SVARs. Journal of Monetary Economics, 117, 798–815.
• Azzalini, A., & Capitanio, A. (1999). Statistical applications of the multivariate skew normal distribution. Journal of the Royal Statistical Society: Series B, 61(3), 579–602.
• Waggoner, D. F., & Zha, T. (1999). Conditional forecasts in dynamic multivariate models. Review of Economics and Statistics, 81(4), 639–651.