cultivars.multivariate.structural.non_gaussian#
Non-Gaussian identification: independence does what zeros used to.
Every impact matrix consistent with the innovation covariance is \(B = L Q\) for the Cholesky factor \(L\) and some orthogonal \(Q\), and the Gaussian distribution is the only one a rotation cannot leave: rotate two independent normal shocks and the result is two independent normal shocks, which is precisely why every Gaussian scheme in this package needs an economic restriction to choose \(Q\). Drop normality and the symmetry breaks. With mutually independent shocks of which at most one is Gaussian, \(Q\) is identified up to column order and sign from the data alone (Comon 1994), and the SVAR literature has built estimators on exactly this fact (Lanne, Meitz and Saikkonen 2017; Gouriéroux, Monfort and Renne 2017). The estimator here is the moment route rather than a parametric likelihood: whiten the innovations, and find the rotation that jointly diagonalizes the third- and fourth-order cumulant slices,
every one of which is diagonal in source coordinates when the sources are independent. Third-order slices carry skewness and fourth-order slices carry tail weight, so a shock identified through either is reachable, and the search is warm-started at the eigenvectors of the kurtosis-weighted covariance – the FOBI solution – on the same joint-diagonalization surface the heteroskedasticity scheme refines.
Two commitments shape the surface. First, the data identify and the user
labels. The recovered shocks are statistical objects, ordered by how
non-Gaussian they are and named shock1 ... shockk unless the user says
otherwise; an economic name for any of them is a claim to be argued from
outside the model, and the restriction note on every summary says so.
Second, the identification condition is checked, not assumed. Each shock’s
skewness and excess kurtosis are reported against their sampling noise, a
count of statistically Gaussian shocks above one is flagged WEAK, and the
residual off-diagonal energy after joint diagonalization is reported as the
data’s verdict on independence itself – the assumption that is strictly
stronger than uncorrelatedness, and the one common volatility across shocks
would violate.
Layout. NonGaussianSVAR is an _IdentificationModel from
_internals; it whitens through _lower_cholesky in _core, builds
the slices with _cumulant_slices in _core, and minimizes the shared
_CoDiagonalObjective from _internals with _solve, then packages
the rotation into SVARResult from
zero_restrictions with
scheme="non-Gaussian". The second-moment route to the same kind of
statistical identification is
heteroskedacity.
References
Comon, P. (1994). Independent component analysis, a new concept? Signal Processing, 36(3), 287-314.
Cardoso, J.-F., & Souloumiac, A. (1993). Blind beamforming for non-Gaussian signals. IEE Proceedings F, 140(6), 362-370.
Lanne, M., Meitz, M., & Saikkonen, P. (2017). Identification and estimation of non-Gaussian structural vector autoregressions. Journal of Econometrics, 196(2), 288-304.
Gouriéroux, C., Monfort, A., & Renne, J.-P. (2017). Statistical inference for independent component analysis: Application to structural VAR models. Journal of Econometrics, 196(1), 111-126.
Example
The same bivariate system identified twice, once with Laplace and uniform shocks and once with Gaussian ones. The first recovers the impact columns from independence alone; the second returns a rotation the data could not have chosen, and says so:
>>> import numpy as np
>>> from cultivars.multivariate.reduced_form.vector_autoregression import VAR
>>> rng = np.random.default_rng(0)
>>> B = np.array([[1.0, 0.5], [0.3, 1.0]])
>>> def simulate(eps):
... y = np.zeros((602, 2))
... for t in range(1, 602):
... y[t] = 0.4 * y[t - 1] + B @ eps[t]
... return VAR(y, order=1).fit()
>>> heavy = rng.laplace(size=602) / np.sqrt(2)
>>> flat = rng.uniform(-np.sqrt(3), np.sqrt(3), size=602)
>>> sharp = NonGaussianSVAR(simulate(np.column_stack([heavy, flat]))).identify()
>>> blunt = NonGaussianSVAR(simulate(rng.standard_normal((602, 2)))).identify()
>>> dict(sharp.diagnostics)["Identification"], dict(blunt.diagnostics)["Identification"]
('at most one Gaussian shock', 'WEAK: 2 shocks statistically Gaussian')
>>> unit = lambda v: v / np.linalg.norm(v)
>>> bool(np.abs(unit(sharp.impact[:, 0]) - unit(B[:, 0])).max() < 0.07)
True
Classes#
Identification from shock independence and non-Gaussianity, Comon (1994). |