from __future__ import annotations
import warnings
from dataclasses import dataclass, field
from typing import Optional, Union
import numpy as np
import numpy.typing as npt
from scipy import linalg as sla
from scipy.stats import chi2
import pandas as pd
def _to_numpy(x: Union[npt.NDArray[np.floating], pd.Series, pd.DataFrame]) -> npt.NDArray[np.floating]:
"""Return ``x`` as a NumPy array with ``float`` dtype.
Args:
x: Array-like input (NumPy, Series, or DataFrame).
Returns:
ndarray: Dense ``float`` array view/copy of ``x``.
"""
return x.to_numpy(dtype=float) if isinstance(x, (pd.Series, pd.DataFrame)) else np.asarray(x, dtype=float)
def _wrap_matrix_like(
arr: npt.NDArray[np.floating], template: Union[npt.NDArray[np.floating], pd.DataFrame]
) -> Union[npt.NDArray[np.floating], pd.DataFrame]:
"""Wrap a matrix to match the pandas/NumPy type of ``template``.
Args:
arr: Matrix to wrap.
template: Reference object that determines the return type.
Returns:
ndarray or pd.DataFrame: ``arr`` with the same container type as ``template``.
"""
if isinstance(template, pd.DataFrame):
return pd.DataFrame(arr, index=template.index, columns=template.columns)
return arr
def _cholesky_pd(
mat: npt.NDArray[np.floating], jitter: float = 1e-12, *, max_attempts: int = 6
) -> npt.NDArray[np.floating]:
r"""Computes a Cholesky decomposition that is robust to numerical errors.
The theoretical framework for robust Bayesian allocation often requires
operations like spectral or Cholesky decompositions of covariance matrices,
which must be positive-definite (PD). However, matrices estimated
from real-world data may fail to be perfectly PD due to floating-point
inaccuracies or estimation errors.
This function attempts to compute the Cholesky decomposition of `mat`. If
`mat` is not PD, it repeatedly adds a multiple of the identity matrix
(`jitter` * I, inflated geometrically) to the matrix and retries. This is a
standard numerical stabilization technique that escalates the diagonal
adjustment until the decomposition succeeds or a maximum number of attempts
is reached.
Args:
mat: The square matrix (N x N) for which to compute the Cholesky
decomposition.
jitter: A small positive constant to seed the diagonal adjustment if
the matrix is not positive-definite. Defaults to 1e-12.
max_attempts: Maximum number of jitter escalations before giving up.
Returns:
The lower Cholesky factor of `mat` (or `mat` + `jitter` * I).
Raises:
ValueError: If `mat` is not a square matrix.
RuntimeError: If the matrix is not positive-definite even after
adding jitter.
"""
mat = np.asarray(mat, dtype=float)
if mat.ndim != 2 or mat.shape[0] != mat.shape[1]:
raise ValueError("Input to _cholesky_pd must be a square matrix.")
if not np.allclose(mat, mat.T, rtol=1e-8, atol=1e-10):
warnings.warn(
"Input matrix is not symmetric (max asymmetry: "
f"{float(np.max(np.abs(mat - mat.T))):.2e}). Symmetrising.",
UserWarning, stacklevel=3,
)
mat = 0.5 * (mat + mat.T)
identity = np.eye(mat.shape[0])
min_jitter = float(jitter) if jitter > 0 else 1e-12
attempts = 0
current_jitter = 0.0
while True:
candidate = mat if current_jitter == 0.0 else mat + current_jitter * identity
try:
return sla.cholesky(candidate, lower=True, check_finite=False)
except sla.LinAlgError as exc:
if attempts >= max_attempts:
raise RuntimeError(
"Matrix not positive-definite after repeated jitter attempts."
) from exc
current_jitter = min_jitter if current_jitter == 0.0 else current_jitter * 10.0
attempts += 1
warnings.warn(
f"Matrix not positive-definite; retrying with jitter={current_jitter:.1e} (attempt {attempts}).",
RuntimeWarning,
)
[docs]
def chi2_quantile(p: float, dof: int, sqrt: bool = False) -> float:
r"""Computes the quantile of the chi-square (\chi^2) distribution.
This function is used to determine the size of the uncertainty ellipsoids
for the market parameters (mean and covariance). The size is
determined by a radius factor, :math:`q`, which is set according to a
quantile of the chi-square distribution.
For the mean vector :math:`\mu`, under the assumption that its posterior
distribution is normal, the squared Mahalanobis distance is chi-square
distributed. The radius factor squared, :math:`q_\mu^2`, is set using a
quantile of the :math:`\chi^2_N` distribution.
.. math::
q_\mu^2 = Q_{\chi_N^2}(p_\mu)
For the covariance matrix :math:`\Sigma`, a similar approach is used based on a
heuristic argument that the Mahalanobis distance behaves like a
:math:`\chi^2` distribution with :math:`N(N+1)/2` degrees of freedom
.
.. math::
q_\Sigma^2 = Q_{\chi_{N(N+1)/2}^2}(p_\Sigma)
Args:
p: The probability level (0 < p < 1) for the quantile.
dof: The degrees of freedom for the chi-square distribution.
sqrt: If True, returns the square root of the quantile. Defaults to False.
Returns:
The chi-square quantile :math:`Q_{\chi^2}(p)` or :math:`\sqrt{Q_{\chi^2}(p)}`
if `sqrt` is True.
Raises:
ValueError: If `p` is not strictly between 0 and 1.
"""
if not (0.0 < p < 1.0):
raise ValueError("Probability `p` must be strictly between 0 and 1.")
q_val = chi2.ppf(p, dof)
return float(np.sqrt(q_val)) if sqrt else float(q_val)
[docs]
@dataclass()
class NIWParams:
r"""A container for the parameters of a Normal-Inverse-Wishart (NIW) posterior distribution. :no-inheritance:
These parameters are the result of a Bayesian update, combining an
NIW prior with market data, as detailed in Meucci (2005). The
formulas for these posterior parameters are given in Eqs. (11)-(14).
"""
T1: int
#: The posterior pseudo-observations for the mean (:math:`T_1`),
#: representing the updated confidence in the mean estimate.
mu1: Union[npt.NDArray[np.floating], "pd.Series[np.floating]"]
#: The posterior mean vector (:math:`\mu_1`), which is the updated
#: estimate of the expected returns.
nu1: int
#: The posterior degrees of freedom for the covariance
#: (:math:`\nu_1`), representing the updated confidence in the
#: covariance estimate.
sigma1: Union[npt.NDArray[np.floating], "pd.DataFrame[np.floating]"]
#: The posterior scale matrix for the covariance (:math:`\Sigma_1`),
#: which is the updated scale matrix of the Inverse-Wishart
#: distribution.
[docs]
class NIWPosterior:
r"""Computes and manages Normal-Inverse-Wishart (NIW) posterior parameters.
This class implements the Bayesian update rules for an NIW distribution,
which is the conjugate prior for a multivariate normal likelihood with
unknown mean and covariance. The methodology follows Section 3
of Meucci (2005). It provides methods to calculate posterior
parameters, classical-equivalent estimators, and factors used in robust
Bayesian asset allocation.
The model assumes that asset returns are independently and identically
distributed according to a normal distribution. The investor's
prior knowledge is modeled as an NIW distribution.
How to Use:
1. Initialize the :class:`NIWPosterior` object with prior parameters:
* ``prior_mu`` (:math:`\mu_0`): The prior estimate for the mean vector.
* ``prior_sigma`` (:math:`\Sigma_0`): The prior scale matrix for the covariance.
* ``t0`` (:math:`T_0`): The confidence in ``prior_mu``, expressed as a
pseudo-count of observations.
* ``nu0`` (:math:`\nu_0`): The confidence in ``prior_sigma``, expressed as a
pseudo-count of observations.
2. Call the :meth:`update` method with sample statistics from observed data:
* ``sample_mu`` (:math:`\hat{\mu}`): The mean vector from the data.
* ``sample_sigma`` (:math:`\hat{\Sigma}`): The covariance matrix from the data.
* ``n_obs`` (:math:`T`): The number of observations in the data sample.
3. The :meth:`update` method returns an :class:`NIWParams` object with the
posterior parameters (:math:`T_1, \mu_1, \nu_1, \Sigma_1`), which are
a blend of the prior and the market data.
4. Use accessor methods like :meth:`get_mu_ce`, :meth:`get_S_mu`, etc., to
retrieve various quantities derived from the posterior distribution.
Attributes:
prior_mu: The prior mean vector (:math:`\mu_0`).
prior_sigma: The prior scale matrix (:math:`\Sigma_0`).
t0: The prior pseudo-count for the mean (:math:`T_0`).
nu0: The prior pseudo-count for the covariance (:math:`\nu_0`).
N: The number of assets.
_asset_index: Stores pandas.Index if pandas objects are used.
_posterior: Stores the computed posterior parameters.
"""
def __init__(
self,
prior_mu: Union[npt.NDArray[np.floating], "pd.Series[np.floating]"],
prior_sigma: Union[npt.NDArray[np.floating], "pd.DataFrame[np.floating]"],
t0: int,
nu0: int,
) -> None:
r"""Initializes the NIWPosterior object with prior parameters.
The prior parameters :math:`(\mu_0, \Sigma_0)` represent the investor's
experience, while :math:`(T_0, \nu_0)` represent their confidence in
that experience.
Args:
prior_mu: 1D array (length N) of prior means (:math:`\mu_0`).
prior_sigma: 2D array (N x N) of the prior scale matrix (:math:`\Sigma_0`).
t0: Prior pseudo-count for the mean (:math:`T_0`). Must be > 0.
nu0: Prior pseudo-count for the covariance (:math:`\nu_0`). Must be >= 0.
Raises:
ValueError: If input parameters have inconsistent shapes or
invalid values.
"""
# Detect pandas inputs
self._pandas = False
self._asset_index: Optional[pd.Index] = None
if isinstance(prior_mu, pd.Series):
self._pandas = True
self._asset_index = prior_mu.index.copy()
self.prior_mu: npt.NDArray[np.floating] = prior_mu.values.astype(float)
else:
self.prior_mu = np.asarray(prior_mu, dtype=float)
if self.prior_mu.ndim != 1:
raise ValueError("`prior_mu` must be a 1D array or pandas.Series.")
self.N: int = self.prior_mu.size
if self.N == 0:
raise ValueError("`prior_mu` cannot be empty; N must be > 0.")
if isinstance(prior_sigma, pd.DataFrame):
if self._asset_index is None:
# If prior_mu was not pandas but prior_sigma is, infer index
self._pandas = True
self._asset_index = prior_sigma.index.copy()
else:
# Ensure indices match
if not prior_sigma.index.equals(self._asset_index) or not prior_sigma.columns.equals(self._asset_index):
raise ValueError("Index/columns of prior_sigma must match index of prior_mu.")
self.prior_sigma: npt.NDArray[np.floating] = prior_sigma.values.astype(float)
else:
self.prior_sigma = np.asarray(prior_sigma, dtype=float)
if self.prior_sigma.ndim != 2 or self.prior_sigma.shape != (self.N, self.N):
raise ValueError(f"`prior_sigma` must be a square matrix of shape ({self.N}, {self.N}), or a pandas.DataFrame with matching index/columns.")
if t0 <= 0:
raise ValueError("`t0` (prior pseudo-count for mean) must be strictly positive.")
if nu0 < 0:
raise ValueError("`nu0` (prior pseudo-count for covariance) must be non-negative.")
_ = _cholesky_pd(self.prior_sigma) # Ensure prior_sigma is PD or near-PD
self.t0: int = int(t0)
self.nu0: int = int(nu0)
self._posterior: Optional[NIWParams] = None
@property
def is_updated(self) -> bool:
"""Whether the posterior has been computed via ``update()``."""
return self._posterior is not None
def _wrap_vector(self, arr: np.ndarray) -> Union[np.ndarray, pd.Series]:
"""Wrap a 1D array as pandas Series if pandas mode is active."""
if self._pandas and self._asset_index is not None:
return pd.Series(arr, index=self._asset_index)
return arr
def _wrap_matrix(self, arr: np.ndarray) -> Union[np.ndarray, pd.DataFrame]:
"""Wrap a 2D array as pandas DataFrame if pandas mode is active."""
if self._pandas and self._asset_index is not None:
return pd.DataFrame(arr, index=self._asset_index, columns=self._asset_index)
return arr
[docs]
def update(
self,
sample_mu: Union[npt.NDArray[np.floating], "pd.Series[np.floating]"],
sample_sigma: Union[npt.NDArray[np.floating], "pd.DataFrame[np.floating]"],
n_obs: int,
) -> NIWParams:
r"""Updates the posterior parameters using sample statistics.
This method implements the Bayesian update rules for the NIW parameters
as given by Eqs. (11)-(14) in Meucci (2005).
.. math::
\begin{align*}
T_1 &= T_0 + T \\
\mu_1 &= \frac{T_0\mu_0 + T\hat{\mu}}{T_1} \\
\nu_1 &= \nu_0 + T \\
\Sigma_1 &= \frac{1}{\nu_1} \left[ \nu_0\Sigma_0 + T\hat{\Sigma} + \frac{(\mu_0 - \hat{\mu})(\mu_0 - \hat{\mu})'}{\frac{1}{T} + \frac{1}{T_0}} \right]
\end{align*}
The resulting posterior parameters blend the investor's prior with
information from the market, with the balance determined by the
relative confidence levels (:math:`T_0, \nu_0`) versus the amount of
data (:math:`T`).
Args:
sample_mu: 1D array (length N) of sample means (:math:`\hat{\mu}`).
sample_sigma: 2D array (N x N) of the sample covariance matrix
(:math:`\hat{\Sigma}`).
n_obs: The number of observations in the sample (:math:`T`).
Returns:
A :class:`NIWParams` instance with the updated posterior parameters.
If pandas objects were used in initialization, returns :math:`\mu_1`
as a Series and :math:`\Sigma_1` as a DataFrame.
Raises:
ValueError: If sample statistics have inconsistent shapes or ``n_obs``
is not a positive integer.
"""
if self._posterior is not None:
raise RuntimeError(
"Posterior already computed. Create a new NIWPosterior for a fresh update."
)
# Convert pandas inputs if present
if isinstance(sample_mu, pd.Series):
smu = sample_mu.values.astype(float)
else:
smu = np.asarray(sample_mu, dtype=float)
if smu.ndim != 1 or smu.shape[0] != self.N:
raise ValueError(f"`sample_mu` must be a 1D array or pandas.Series of length {self.N}.")
if n_obs <= 0:
raise ValueError("`n_obs` (number of observations) must be strictly positive.")
if isinstance(sample_sigma, pd.DataFrame):
ssigma = sample_sigma.values.astype(float)
else:
ssigma = np.asarray(sample_sigma, dtype=float)
if ssigma.ndim != 2 or ssigma.shape != (self.N, self.N):
raise ValueError(f"`sample_sigma` must be a square matrix of shape ({self.N}, {self.N}), or pandas.DataFrame with matching index/columns.")
_ = _cholesky_pd(ssigma) # Ensure sample_sigma is PD or near-PD
# Compute posterior scalars
T1 = self.t0 + n_obs
nu1 = self.nu0 + n_obs
mu1_array = (self.t0 * self.prior_mu + n_obs * smu) / T1
cross_term_weight_denominator = (1.0 / n_obs + 1.0 / self.t0)
diff_mu = self.prior_mu - smu
outer_prod_diff_mu = np.outer(diff_mu, diff_mu)
cross_term_weighted = outer_prod_diff_mu / cross_term_weight_denominator
sigma1_numerator = (self.nu0 * self.prior_sigma
+ n_obs * ssigma
+ cross_term_weighted)
if nu1 <= 0:
raise ValueError("Posterior degrees of freedom \nu_1 must be positive.")
sigma1_array = sigma1_numerator / nu1
_ = _cholesky_pd(sigma1_array) # Ensure \Sigma_1 is PD or near-PD
# Wrap into pandas if appropriate
mu1 = self._wrap_vector(mu1_array)
sigma1 = self._wrap_matrix(sigma1_array)
self._posterior = NIWParams(T1=T1, mu1=mu1, nu1=nu1, sigma1=sigma1)
return self._posterior
[docs]
def get_mu_ce(self) -> Union[npt.NDArray[np.floating], "pd.Series[np.floating]"]:
r"""Computes the classical-equivalent estimator for the mean, :math:`\hat{\mu}_{ce}`.
For the NIW model, this estimator is the posterior mean :math:`\mu_1`,
as defined in Eq. (15).
.. math::
\hat{\mu}_{ce} = \mu_1
Returns:
The posterior mean vector :math:`\mu_1` as a NumPy array or pandas
Series.
Raises:
RuntimeError: If posterior parameters have not been computed via
:meth:`update`.
"""
if self._posterior is None:
raise RuntimeError("Posterior parameters not computed. Call `update()` first.")
return self._wrap_vector(np.asarray(self._posterior.mu1, dtype=float).ravel())
[docs]
def get_S_mu(self) -> Union[npt.NDArray[np.floating], "pd.DataFrame[np.floating]"]:
r"""Computes the scatter matrix :math:`S_{\mu}` for the posterior of :math:`\mu`.
This matrix describes the dispersion of the marginal posterior
distribution of :math:`\mu` and is used to define the location-dispersion
ellipsoid for robust optimization. It is defined in Eq. (16):
.. math::
S_{\mu} = \frac{1}{T_1} \frac{\nu_1}{\nu_1 - 2} \Sigma_1
This computation requires the posterior degrees of freedom :math:`\nu_1 > 2`.
Returns:
The scatter matrix :math:`S_{\mu}` as a NumPy array or pandas DataFrame.
Raises:
RuntimeError: If posterior parameters have not been computed.
ValueError: If :math:`\nu_1 \le 2`, as the scatter matrix is not
defined.
"""
if self._posterior is None:
raise RuntimeError("Posterior parameters not computed. Call `update()` first.")
if self._posterior.nu1 <= 2:
raise ValueError(
f"Posterior nu1={self._posterior.nu1} must be > 2 for S_mu "
f"(nu1 = nu0 + T = {self._posterior.nu1}). "
"Increase prior nu0 or provide more observations."
)
factor = self._posterior.nu1 / (self._posterior.T1 * (self._posterior.nu1 - 2.0))
sigma1_array = _to_numpy(self._posterior.sigma1)
S_mu_array = factor * sigma1_array
return self._wrap_matrix(S_mu_array)
[docs]
def get_sigma_ce(self) -> Union[npt.NDArray[np.floating], "pd.DataFrame[np.floating]"]:
r"""Computes the classical-equivalent estimator for the covariance, :math:`\hat{\Sigma}_{ce}`.
This estimator is a shrunk version of the posterior scale matrix
:math:`\Sigma_1`, as defined in Eq. (17). It serves as the center of
the uncertainty ellipsoid for :math:`\Sigma`.
.. math::
\hat{\Sigma}_{ce} = \frac{\nu_1}{\nu_1 + N + 1} \Sigma_1
Returns:
The classical-equivalent estimator :math:`\hat{\Sigma}_{ce}` as a NumPy
array or pandas DataFrame.
Raises:
RuntimeError: If posterior parameters have not been computed.
ValueError: If :math:`\nu_1 + N + 1 = 0`, which is highly unlikely
with valid inputs.
"""
if self._posterior is None:
raise RuntimeError("Posterior parameters not computed. Call `update()` first.")
denom = self._posterior.nu1 + self.N + 1.0
if denom == 0:
raise ValueError(r"Denominator (\nu_1 + N + 1) for \Sigma_ce is zero.")
factor = self._posterior.nu1 / denom
sigma1_array = _to_numpy(self._posterior.sigma1)
sigma_ce_array = factor * sigma1_array
return self._wrap_matrix(sigma_ce_array)
[docs]
def cred_radius_mu(self, p_mu: float) -> float:
r"""Computes the credibility factor :math:`\gamma_\mu` for the mean's uncertainty.
This factor, :math:`\gamma_\mu`, appears in the simplified robust
mean-variance optimization problem (Eq. 19). It scales the
portfolio's posterior standard deviation to penalize for estimation
risk in the mean vector. Its formula is given by Eq. (20):
.. math::
\gamma_\mu = \sqrt{ \frac{q_\mu^2}{T_1} \frac{\nu_1}{\nu_1 - 2} }
where :math:`q_\mu^2 = Q_{\chi^2_N}(p_{mu})` is the squared radius factor
from the chi-square distribution.
Args:
p_mu: The confidence level for :math:`\mu` (0 < p_mu < 1), which
reflects aversion to estimation risk.
Returns:
The credibility factor :math:`\gamma_\mu`.
Raises:
RuntimeError: If posterior parameters have not been computed.
ValueError: If :math:`\nu_1 \le 2` or `p_mu` is not in (0,1).
"""
if self._posterior is None:
raise RuntimeError("Posterior parameters not computed. Call `update()` first.")
if self._posterior.nu1 <= 2:
raise ValueError(r"Posterior \nu_1 must be greater than 2 for \gamma_\mu calculation.")
if not (0.0 < p_mu < 1.0):
raise ValueError("Confidence level p_mu must be between 0 and 1 (exclusive).")
q_mu_squared = chi2_quantile(p_mu, self.N, sqrt=False)
term_T1 = self._posterior.T1
term_nu1 = self._posterior.nu1
gamma_mu = np.sqrt((q_mu_squared / term_T1) * (term_nu1 / (term_nu1 - 2.0)))
return gamma_mu
[docs]
def cred_radius_sigma_factor(self, p_sigma: float) -> float:
r"""Computes the scaling factor for the worst-case portfolio variance.
In the robust framework, the maximum possible variance of a portfolio
within the uncertainty ellipsoid for :math:`\Sigma` is not simply
:math:`w'\Sigma_1 w`, but a scaled version of it. This method computes
that scaling factor, which we can call :math:`C_\Sigma`. The derivation is
shown in Appendix 7.2, and the final result is presented in the
maximization step in Eq. (47):
.. math::
\max_{\Sigma \in \Theta_\Sigma} w'\Sigma w = \underbrace{ \left[ \frac{\nu_1}{\nu_1 + N + 1} + \sqrt{\frac{2\nu_1^2 q_\Sigma^2}{(\nu_1 + N + 1)^3}} \right] }_{C_\Sigma} (w'\Sigma_1 w)
where :math:`q_\Sigma^2 = Q_{\chi^2_{dof}}(p_\Sigma)` with
:math:`dof = N(N+1)/2`.
Args:
p_sigma: The confidence level for :math:`\Sigma` (0 < p_sigma < 1),
reflecting aversion to estimation risk.
Returns:
The credibility factor :math:`C_\Sigma` for scaling the portfolio variance.
Raises:
RuntimeError: If posterior parameters have not been computed.
ValueError: If `p_sigma` is not in (0,1) or internal terms are invalid.
"""
if self._posterior is None:
raise RuntimeError("Posterior parameters not computed. Call `update()` first.")
if not (0.0 < p_sigma < 1.0):
raise ValueError("Confidence level p_sigma must be between 0 and 1 (exclusive).")
nu1 = self._posterior.nu1
denom_cubed_base = nu1 + self.N + 1.0
if denom_cubed_base <= 0:
raise ValueError(r"Term (\nu_1 + N + 1) must be positive for C_\Sigma calculation.")
dof = self.N * (self.N + 1) // 2
q_sigma_squared_val = chi2_quantile(p_sigma, dof, sqrt=False)
term1 = nu1 / denom_cubed_base
numerator_term2 = 2.0 * nu1**2 * q_sigma_squared_val
denominator_term2 = denom_cubed_base**3
if denominator_term2 == 0:
raise ValueError(r"Denominator for sqrt term in C_\Sigma is zero.")
term2_arg = numerator_term2 / denominator_term2
if term2_arg < 0:
warnings.warn(
f"Numerical noise in C_Sigma: sqrt argument is {term2_arg:.2e}; clamping to zero.",
RuntimeWarning, stacklevel=2,
)
term2_arg = 0.0
term2 = np.sqrt(term2_arg)
C_sigma = term1 + term2
return C_sigma
[docs]
@dataclass(frozen=True)
class RobustBayesPosterior:
"""Bundle NIW posterior moments with mean-uncertainty helpers.
The Normal-Inverse-Wishart (NIW) posterior implies a covariance of the mean
:math:`S_\\mu` given by :cite:p:`meucci2005robust`:
.. math::
S_\\mu = \\frac{\\nu_1}{T_1 (\\nu_1 - 2)} \\Sigma_1.
This object exposes :math:`S_\\mu` and convenience transforms into
horizon-scaled log- and simple-return units, which are required by robust
optimisation routines that penalise mean uncertainty.
Key fields:
mu: Posterior mean vector (same units as inputs).
sigma: Posterior covariance matrix (same units as inputs).
s_mu: Posterior covariance of the mean (:math:`S_\\mu`).
"""
mu: Union[npt.NDArray[np.floating], pd.Series]
sigma: Union[npt.NDArray[np.floating], pd.DataFrame]
s_mu: Union[npt.NDArray[np.floating], pd.DataFrame]
_nu1: Optional[float] = field(default=None, repr=False)
_n_assets: Optional[int] = field(default=None, repr=False)
_T1: Optional[float] = field(default=None, repr=False)
[docs]
@classmethod
def from_niw(
cls,
*,
prior_mu: Union[npt.NDArray[np.floating], pd.Series],
prior_sigma: Union[npt.NDArray[np.floating], pd.DataFrame],
t0: int,
nu0: int,
sample_mu: Union[npt.NDArray[np.floating], pd.Series],
sample_sigma: Union[npt.NDArray[np.floating], pd.DataFrame],
n_obs: int,
) -> "RobustBayesPosterior":
"""Construct the posterior bundle from NIW inputs.
Args:
prior_mu: Prior mean vector.
prior_sigma: Prior covariance matrix.
t0: Prior strength for the mean.
nu0: Prior degrees of freedom for the covariance.
sample_mu: Sample mean vector.
sample_sigma: Sample covariance matrix.
n_obs: Number of observations.
Returns:
RobustBayesPosterior: Posterior bundle with mean uncertainty.
"""
niw = NIWPosterior(
prior_mu=prior_mu,
prior_sigma=prior_sigma,
t0=int(t0),
nu0=int(nu0),
)
niw.update(sample_mu=sample_mu, sample_sigma=sample_sigma, n_obs=int(n_obs))
return cls(
mu=niw.get_mu_ce(),
sigma=niw.get_sigma_ce(),
s_mu=niw.get_S_mu(),
_nu1=float(niw._posterior.nu1),
_n_assets=int(niw.N),
_T1=float(niw._posterior.T1),
)
[docs]
def mean_uncertainty_cov_log(
self,
*,
annualization_factor: float = 1.0,
) -> Union[npt.NDArray[np.floating], pd.DataFrame]:
"""Return mean-uncertainty covariance in log-return units.
The mean uncertainty scales with the square of the horizon:
.. math::
S_{\\mu,\\,h} = h^2\\, S_\\mu.
"""
s_mu_np = _to_numpy(self.s_mu)
scaled = s_mu_np * (annualization_factor ** 2)
return _wrap_matrix_like(scaled, self.s_mu)
@property
def sigma1(self) -> Union[npt.NDArray[np.floating], pd.DataFrame]:
"""Return the raw posterior scale matrix Sigma_1 (Meucci Eq. 7.31).
``sigma`` stores the classical-equivalent Sigma_ce = nu1/(nu1+N+1) * Sigma_1.
This property recovers Sigma_1 = (nu1+N+1)/nu1 * Sigma_ce for use
in the robust Bayesian formulation (Meucci Eq. 9.155).
Returns:
Posterior scale matrix Sigma_1.
Raises:
ValueError: If the posterior was not constructed via ``from_niw``.
"""
if self._nu1 is None or self._n_assets is None:
raise ValueError("sigma1 requires NIW posterior params; use RobustBayesPosterior.from_niw().")
factor = (self._nu1 + self._n_assets + 1.0) / self._nu1
sigma_ce_np = _to_numpy(self.sigma)
sigma1_np = factor * sigma_ce_np
return _wrap_matrix_like(sigma1_np, self.sigma)
[docs]
def cred_radius_mu(self, p_mu: float = 0.75) -> float:
"""Credibility radius for the mean uncertainty ellipsoid (Meucci Eq. 9.156).
Args:
p_mu: Probability level for the chi-square quantile. Defaults to 0.75.
Returns:
float: Credibility factor gamma_mu.
Raises:
ValueError: If NIW parameters are not available.
"""
if self._nu1 is None or self._n_assets is None or self._T1 is None:
raise ValueError("cred_radius_mu requires NIW params; use from_niw().")
q_mu_sq = chi2_quantile(p_mu, self._n_assets, sqrt=False)
return float(np.sqrt((q_mu_sq / self._T1) * (self._nu1 / (self._nu1 - 2.0))))
[docs]
def cred_radius_sigma_factor(self, p_sigma: float = 0.75) -> float:
"""Scaling factor for worst-case variance under Sigma uncertainty (Meucci Eq. 9.157).
Args:
p_sigma: Probability level for the chi-square quantile. Defaults to 0.75.
Returns:
float: Scaling factor C_sigma.
Raises:
ValueError: If NIW parameters are not available.
"""
if self._nu1 is None or self._n_assets is None:
raise ValueError("cred_radius_sigma_factor requires NIW params; use from_niw().")
denom = self._nu1 + self._n_assets + 1.0
dof = self._n_assets * (self._n_assets + 1) // 2
q_sig_sq = chi2_quantile(p_sigma, dof, sqrt=False)
term1 = self._nu1 / denom
term2_arg = 2.0 * self._nu1**2 * q_sig_sq / denom**3
if term2_arg < 0:
warnings.warn(
f"Numerical noise in C_Sigma: clamping sqrt arg {term2_arg:.2e} to zero.",
RuntimeWarning, stacklevel=2,
)
term2_arg = 0.0
return float(term1 + np.sqrt(term2_arg))
[docs]
def mean_uncertainty_cov_simple(
self,
*,
annualization_factor: float = 1.0,
exact: bool = False,
) -> Union[npt.NDArray[np.floating], pd.DataFrame]:
"""Return mean-uncertainty covariance in simple-return units.
Assumes ``mu`` and ``s_mu`` are expressed in log-return units. The method
annualises and applies a delta approximation to map to simple returns:
.. math::
r \\approx \\exp(\\mu_g) - 1, \\quad
S_{\\mu,r} \\approx J\\,S_{\\mu,g}\\,J^{\\top}, \\quad
J = \\operatorname{diag}(\\exp(\\mu_g)).
Args:
annualization_factor: Horizon scaling factor (e.g. 12 for monthly -> annual).
exact: When ``True``, use the exact delta-method Jacobian
``J = diag(exp(mu_g + 0.5 * diag(Sigma_g)))`` which accounts for
the Jensen inequality correction. Default ``False`` uses the
simpler ``J = diag(exp(mu_g))``.
"""
mu_np = _to_numpy(self.mu)
s_mu_np = _to_numpy(self.s_mu)
sigma_np = _to_numpy(self.sigma)
mu_log = mu_np * annualization_factor
s_mu_log = s_mu_np * (annualization_factor ** 2)
if exact:
sigma_log = sigma_np * (annualization_factor ** 2)
scale = np.exp(mu_log + 0.5 * np.diag(sigma_log))
else:
scale = np.exp(mu_log)
cov_simple = (scale[:, None] * s_mu_log) * scale[None, :]
return _wrap_matrix_like(cov_simple, self.s_mu)