bayesreconpy.linear.conditioning#

bayesreconpy.linear.gaussian#

Gaussian linear reconciliation.

bayesreconpy.linear.gaussian.reconc_gaussian(A, base_forecasts_mu, base_forecasts_Sigma)#

Analytical reconciliation of Gaussian base forecasts. Reconciles forecasts using a closed-form computation for Gaussian base forecasts.

Parameters:

  • A (ndarray) – Aggregation matrix with shape (n_upper, n_bottom). Each column represents a bottom-level forecast, and each row represents an upper-level forecast. A value of 1 in A[i, j] indicates that the bottom-level forecast j contributes to the upper-level forecast i.

  • base_forecasts_mu (List[float]) –

    A 1D list containing the means of the base forecasts. The order is:

    • First, the upper-level means (in the order of rows in A).

    • Then, the bottom-level means (in the order of columns in A).

  • base_forecasts_Sigma (ndarray) – A 2D covariance matrix representing the uncertainties in the base forecasts. The order of rows and columns must match the order of base_forecasts_mu. Shape: (n, n), where n = n_upper + n_bottom.

Returns:

A dictionary containing:

  • ’bottom_reconciled_mean’numpy.ndarray

    A 1D array of the reconciled means for the bottom-level forecasts. Shape: (n_bottom,).

  • ’bottom_reconciled_covariance’numpy.ndarray

    A 2D covariance matrix of the reconciled bottom-level forecasts. Shape: (n_bottom, n_bottom).

Return type:

Dict[str, ndarray]

Notes

  • The function assumes that base forecasts follow a multivariate Gaussian distribution.

  • The covariance matrix base_forecasts_Sigma should be symmetric and positive semi-definite.

  • The order of elements in base_forecasts_mu and rows/columns in base_forecasts_Sigma is critical: first the upper-level forecasts (in the order of rows in A), followed by the bottom-level forecasts (in the order of columns in A).

  • The function returns only the reconciled parameters for the bottom-level forecasts. Reconciled parameters for upper-level forecasts and the entire hierarchy can be derived using the reconciliation matrix A.

Examples

Example 1: Minimal hierarchy with Gaussian base forecasts
>>> A = np.array([
...     [1, 0, 0],
...     [0, 1, 1]
... ])
>>> base_forecasts_mu = [9.0, 2.0, 4.0]
>>> base_forecasts_Sigma = np.diag([9.0, 4.0, 4.0])
>>> result = reconc_gaussian(A, base_forecasts_mu, base_forecasts_Sigma)
>>> print(result['bottom_reconciled_mean'])
[2.5, 4.0]
>>> print(result['bottom_reconciled_covariance'])
[[2.25, 1.5 ],
 [1.5 , 2.0 ]]

References

  • Corani, G., Azzimonti, D., Augusto, J.P.S.C., Zaffalon, M. (2021). Probabilistic Reconciliation of Hierarchical Forecast via Bayes’ Rule. ECML PKDD 2020. Lecture Notes in Computer Science, vol 12459. https://doi.org/10.1007/978-3-030-67664-3_13

  • Zambon, L., Agosto, A., Giudici, P., Corani, G. (2024). Properties of the reconciled distributions for Gaussian and count forecasts. International Journal of Forecasting (in press). https://doi.org/10.1016/j.ijforecast.2023.12.004

bayesreconpy.linear.projection#

Projection-based linear reconciliation methods.

bayesreconpy.linear.projection.reconc_mint(A, base_forecasts, res=None, W=None, samples=False)#

Reconciles base forecasts using the MinT (Minimum Trace) approach with a shrinkage-based covariance estimator.

The MinT method adjusts base forecasts to ensure coherence with a given aggregation structure, minimizing the total variance of the reconciliation errors using the residual covariance matrix.

Parameters:#

Anp.ndarray

The aggregation constraint matrix (typically of shape [n_agg_levels, n_total_series]), defining how the bottom-level time series aggregate into upper levels.

base_forecastsnp.ndarray

The base forecasts to be reconciled: - of shape [n_total_series, n_time, n_samples] if samples=True, or - of shape [n_total_series, n_time] if samples=False.

resnp.ndarray

Residuals (forecast errors) from a previous model, used to estimate the covariance matrix. Should be of shape [n_total_series, n_time].

Wnp.ndarray, optional

Pre-computed covariance/weight matrix (W_h). If provided, ‘res’ is ignored.

samplesbool, optional (default=False)

Indicates whether the base forecasts include multiple samples. If True, reconciliation is applied sample-by-sample along the third axis.

Returns:#

:

y_tilde_meannp.ndarray

The reconciled forecast means, with the same shape as base_forecasts.

y_tilde_varnp.ndarray or list of np.ndarray

The reconciled forecast variances: - A single matrix of shape [n_total_series, n_total_series] if samples=False. - A list of such matrices, one per sample, if samples=True.

bayesreconpy.linear.projection.reconc_ols(A, base_forecasts, samples=False)#

Reconciles base forecasts using Ordinary Least Squares (OLS) based on a summing matrix.

This function applies forecast reconciliation using the OLS method, ensuring that the reconciled forecasts are coherent with the aggregation structure defined by the matrix A.

Parameters:#

Anp.ndarray

The aggregation constraint matrix (usually of shape [n_agg_levels, n_total_series]). It defines how bottom-level time series are aggregated into upper levels.

base_forecastsnp.ndarray

The base forecasts before reconciliation. Should be: - of shape [n_total_series, n_time, n_samples] if samples=True, or - of shape [n_total_series, n_time] if samples=False.

samplesbool, optional (default=False)

Indicates whether the base forecasts contain multiple samples (i.e., are probabilistic). If True, reconciliation is applied to each sample individually along the third axis.

Returns:#

:

y_tilde_meannp.ndarray

The reconciled forecasts: - of shape [n_total_series, n_time, n_samples] if samples=True, or - of shape [n_total_series, n_time] if samples=False.

bayesreconpy.nonlinear.conditioning#

bayesreconpy.nonlinear.conditioning.reconc_nl_is(f_u, num_samples=10000, n_free=None, seed=None, joint_mean=None, joint_cov=None, eps_cov=1e-09)#

Nonlinear Probabilistic Reconciliation using Importance Sampling (IS).

Uses importance sampling to reconcile forecasts by conditioning on a nonlinear manifold. The method draws samples from the joint distribution of free and constrained variables, computes importance weights based on the constraint function, and resamples to obtain the reconciled distribution.

Parameters:

  • f_u (Callable[[ndarray], ndarray]) –

    The FTC function mapping free variables to constrained variables.

    • Input: array of shape (N, n_free) - N samples of the n_free free variables.

    • Output: array of shape (n_constrained, N) - constrained variables for each sample.

    The function defines the relationship: constrained = f_u(free).

  • num_samples (int, optional) – Number of samples to draw from the reconciled distribution. Default is 10,000.

  • n_free (int) – Number of free (bottom-level) variables. Required parameter. The joint mean and covariance have dimension n_constrained + n_free.

  • seed (int or None, optional) – Random seed for reproducibility. Default is None.

  • joint_mean (numpy.ndarray) –

    Mean vector of the joint distribution over [constrained, free] variables. Shape: (n_constrained + n_free,).

    The first n_constrained elements are the means of the constrained variables, and the remaining n_free elements are the means of the free variables.

  • joint_cov (numpy.ndarray) – Covariance matrix of the joint distribution over [constrained, free] variables. Shape: (n_constrained + n_free, n_constrained + n_free). Must be symmetric and positive semi-definite.

  • eps_cov (float) – Small regularization value added to the diagonal of covariance matrices to ensure numerical stability. Default is 1e-9.

Returns:

A dictionary containing the reconciled forecasts:

  • ’free_reconciled_samples’: numpy.ndarray

    Reconciled samples for the free variables. Shape: (n_free, num_samples).

  • ’constrained_reconciled_samples’: numpy.ndarray

    Reconciled samples for the constrained variables. Shape: (n_constrained, num_samples).

  • ’reconciled_samples’: numpy.ndarray

    Concatenated reconciled samples for all variables [constrained, free]. Shape: (n_constrained + n_free, num_samples).

Return type:

dict

Raises:

  • ValueError – If joint_mean or joint_cov is None. Both parameters are required.

  • TypeError – If joint_mean is None and accessed before validation.

Notes

  • The function computes weights as: w_j ∝ p(U_pred_j, B_j) / p(B_j) where U_pred_j = f_u(B_j) are the predicted constrained variables.

  • Weights are normalized and used to resample particles via multinomial resampling.

  • The joint distribution should reflect the prior beliefs about both free and constrained variables. Ideally, it captures correlations between hierarchical levels.

  • This method is particularly useful for complex nonlinear relationships where the constraint function cannot be inverted analytically.

Examples

Example: Multiple constraints
>>> n_free = 3
>>> n_constrained = 2
>>> n_total = n_free + n_constrained
>>>
>>> joint_mean = np.ones(n_total) * 2.0
>>> joint_cov = np.eye(n_total) * 0.3
>>>
>>> def f_u(B):
...     # B shape: (N, n_free)
...     # Two constraints: sum of first two frees, and third free
...     return np.vstack([
...         np.sum(B[:, :2], axis=1),    # (N,)
...         B[:, 2]                       # (N,)
...     ])  # Output: (n_constrained, N)
>>>
>>> result = reconc_nl_is(f_u, num_samples=5000, n_free=n_free,
...                        joint_mean=joint_mean, joint_cov=joint_cov, seed=42)
>>> print("Reconciled shape:", result["reconciled_samples"].shape)
Reconciled shape: (5, 5000)
bayesreconpy.nonlinear.conditioning.reconc_nl_ukf(free_base_forecasts, in_type, distr, f_u, constrained_base_forecasts, R, num_samples=10000, seed=None)#

Nonlinear Probabilistic Reconciliation using Unscented Kalman Filter (UKF) Conditioning.

Uses the Unscented Kalman Filter to condition the free forecast distributions on the constrained base forecasts via a user-defined nonlinear constraint function. The method performs probabilistic reconciliation by projecting the joint distribution onto the nonlinear manifold defined by f_u(free) = constrained_base_forecasts.

Parameters:

  • free_base_forecasts (List[Dict[str, Union[float, ndarray]]]) –

    A list of n_free dictionaries, one for each free-level forecast.

    • If in_type[i] == “samples”, free_base_forecasts[i] must contain: * ‘samples’: numpy.ndarray of shape (n_samples,) containing samples from the forecast. * ‘residuals’: numpy.ndarray of shape (n_samples,) containing residuals (samples - mean).

    • If in_type[i] == “params”, free_base_forecasts[i] must be a dictionary containing parameters for the specified distribution in distr[i]: * ‘gaussian’: {“mean”: float, “sd”: float}

  • in_type (List[str]) –

    Specifies the input type for each free forecast. Must have length n_free. Each element is one of:

    • ’samples’: The forecast is provided as samples.

    • ’params’: The forecast is provided as parameters.

    Important: All elements must be the same (either all ‘samples’ or all ‘params’). Mixed input types are not supported and will raise NotImplementedError.

  • distr (List[str]) –

    Specifies the distribution type for each free forecast. Must have length n_free.

    • If in_type[i] == “samples”: distr[i] should be descriptive (e.g., ‘gaussian’, ‘continuous’).

    • If in_type[i] == “params”: distr[i] must be one of: * ‘gaussian’: Gaussian distribution

  • f_u (Callable[[ndarray], ndarray]) –

    The FTC function. Maps free variables to constrained variables via the relationship: constrained = f_u(free).

    The function should: - Accept an array of shape (n_free,) representing the free variables. - Return an array of shape (n_constrained,) representing the constrained variables. - Be differentiable (required by the unscented transform).

  • constrained_base_forecasts (ndarray) – Base forecast for the constrained (upper-level) variables. Shape: (n_constrained,). These are the target values that the reconciled free forecasts should aggregate to (according to the constraint function f_u).

  • R (ndarray) – Measurement noise covariance matrix representing uncertainty in the constrained forecasts. Shape: (n_constrained, n_constrained). Must be symmetric and positive semi-definite.

  • num_samples (int) – Number of samples to draw from the reconciled distribution. Default is 10,000.

  • seed (Optional[int]) – Random seed for reproducibility. Default is None.

Returns:

A dictionary containing the reconciled forecasts:

  • ’free_reconciled_samples’: numpy.ndarray

    Reconciled samples for the free (bottom-level) variables. Shape: (n_free, num_samples). Each column is a sample from the reconciled distribution conditioned on the constrained base forecasts.

Return type:

Dict[str, ndarray]

Raises:

  • NotImplementedError – If in_type contains mixed ‘params’ and ‘samples’ entries.

  • ValueError – If required keys are missing from parameter dictionaries or if sample arrays have incompatible shapes.

  • KeyError – If a required parameter key is missing (e.g., ‘sd’ for Gaussian distribution).

Notes

  • The unscented Kalman filter uses sigma points to approximate the mean and covariance of the nonlinearly transformed distributions, avoiding the need for linearization.

  • All free forecasts must use the same input type. For mixed input types, first convert or reconstruct one type before calling this function.

  • The constraint function f_u should be well-defined for the range of the free forecasts.

  • The measurement noise covariance R should reflect the uncertainty in the constrained base forecasts. Smaller values increase the weight of the constrained forecasts.

Examples

ExampleMultiplicative constraint with samples
>>> rng = np.random.default_rng(42)
>>> n_samples = 500
>>> s0 = rng.normal(2.0, 0.2, n_samples)
>>> s1 = rng.normal(3.0, 0.2, n_samples)
>>>
>>> free_base_forecasts = [
...     {"samples": s0, "residuals": s0 - s0.mean()},
...     {"samples": s1, "residuals": s1 - s1.mean()},
... ]
>>> in_type = ["samples"] * 2
>>> distr = ["gaussian"] * 2
>>>
>>> def multiplicative_constraint(free):
...     # free[0] * free[1] should equal constrained value
...     return free[0] * free[1]
>>>
>>> constrained_base_forecasts = np.array([6.0])
>>> R = 0.1 * np.eye(1)
>>>
>>> result = reconc_nl_ukf(free_base_forecasts, in_type, distr,
...                         multiplicative_constraint, constrained_base_forecasts, R,
...                         num_samples=5000, seed=42)
>>> print("Reconciled shape:", result["free_reconciled_samples"].shape)
Reconciled shape: (2, 5000)

bayesreconpy.nonlinear.projection#

bayesreconpy.nonlinear.projection.reconc_nl_ols(samples, f, W=None, n_iter=10, seed=None)#

Nonlinear Probabilistic Reconciliation via Constrained Projection (OLS).

Projects joint forecast samples onto a nonlinear constraint manifold using iterative constrained optimization. The method minimizes the weighted Euclidean distance from original samples to the constraint manifold f(Z) = 0, ensuring all reconciled samples satisfy the nonlinear constraints exactly.

Parameters:

  • samples (ndarray) – Joint forecast samples to be reconciled. Shape: (n_samples, d) where d is the total number of variables (constrained + free variables combined). Each row represents one sample from the joint distribution.

  • f (Callable[[Array], Array]) –

    Constraint function defining the nonlinear reconciliation manifold.

    • Input: JAX array of shape (d,) representing one sample point.

    • Output: JAX array of shape (k,) where k = n_constraints is the number of constraints. The function defines the constraints as f(Z) = 0.

    The constraints express relationships between variables that reconciled samples must satisfy exactly. For example, hierarchical aggregation constraints: top = middle * b2 becomes top - middle * b2 = 0.

  • W (Optional[ndarray]) –

    Weight matrix for the optimization objective. Shape: (d, d). Default is identity matrix np.eye(d).

    Controls the relative importance of different variables when projecting onto the constraint manifold. Can be used to emphasize certain variables or account for different units of measurement.

    The optimization minimizes: ||W @ (Z_proj - Z_orig)||_2^2 subject to f(Z_proj) = 0.

  • n_iter (int) – Number of Newton iterations for solving the constrained optimization problem. Default is 10. Increase for tighter constraint satisfaction or complex manifolds.

  • seed (Optional[int]) – Random seed for reproducibility. Default is None.

Returns:

A dictionary containing the reconciled forecasts and constraint residuals:

  • ’reconciled_samples’: np.ndarray

    Projected samples satisfying the constraint manifold. Shape: (d, n_samples). Each column is a reconciled sample with f(Z_proj) ≈ 0.

  • ’residuals’: np.ndarray

    Constraint residuals f(Z_proj) for each reconciled sample. Shape: (n_constraints, n_samples). Close to zero for successful reconciliation. Larger residuals indicate that the constraint manifold is difficult to satisfy or the optimization did not converge well.

Return type:

Dict[str, ndarray]

Notes

  • Constraint Satisfaction: Unlike sampling-based methods, reconciled samples exactly (or approximately) satisfy the nonlinear constraints f(Z) = 0. This ensures perfect coherency.

  • Optimization Method: Uses constrained Newton iterations (via JNLR library’s make_solver) to project samples onto the manifold. The method solves a constrained least-squares problem at each iteration.

  • Nonlinearity Support: Handles arbitrary smooth nonlinear constraints (multiplicative, exponential, etc.) without linearization assumptions.

  • Weight Matrix: Use W to handle different variable scales or to weight constraints differently. For example, W = np.diag([1.0, 0.1, 0.1, 0.1, 0.1]) prioritizes getting the first variable correct.

  • Convergence: The n_iter parameter controls the Newton iterations. The method should converge for well-posed problems within 10-20 iterations. Check the residuals to assess convergence.

Raises:

  • ValueError – If constraint function input/output shapes are incompatible with samples.

  • TypeError – If samples are not numpy arrays or constraint function is not callable.

Examples

ExampleMultiple constraints with identity weight matrix
>>> rng = np.random.default_rng(42)
>>> n_samples = 2000
>>>
>>> # Create samples
>>> samples = rng.normal(0, 1, (n_samples, 4))
>>>
>>> # Define two constraints
>>> def f(z):
...     return jnp.array([
...         z[0] + z[1] - z[2],         # constraint 1: z0 + z1 = z2
...         z[1] * z[3] - z[0]          # constraint 2: z1 * z3 = z0
...     ])
>>>
>>> result = reconc_nl_ols(samples, f, n_iter=20, seed=42)
>>> Z_proj = result["reconciled_samples"]
>>> residuals = result["residuals"]
>>>
>>> print("Residuals shape:", residuals.shape)
Residuals shape: (2, 2000)
>>> print("Mean abs residuals per constraint:", np.mean(np.abs(residuals), axis=1))
Mean abs residuals per constraint: [1.2e-04 5.8e-05]

bayesreconpy.core#

Core foundational utilities shared by reconciliation methods.

bayesreconpy.utils#

Split utility modules for validation, sampling, constants, and stats.