Bootstrap inference for the latent VAR of an NMF-VAR model

This function is experimental. The interface may change in future versions: argument names, defaults and the contents of the returned object are not yet stable.

nmfkc.ar.latent.inference attaches a wild bootstrap to the quantities reported by nmfkc.ar.latent: standard errors and intervals for the spectral radius, the long-run mean, the long-run composition and the entries of the latent transition matrices $G_d$.

Usage

nmfkc.ar.latent.inference(object, Y, A, ...)

Arguments

object

An object from nmfkc.ar.latent.

Y, A

The observation and covariate matrices the model was fitted on, i.e. the Y and A returned by nmfkc.ar.

...

Additional arguments:

wild.B: Number of bootstrap replicates. Default 500.
wild.dist: Multiplier distribution, "rademacher" (default) or "exp" (mean-centred).
wild.unit: "column" (default, one multiplier per time point) or "element" (one per cell).
wild.level: Confidence level. Default 0.95.
wild.seed: Seed for the bootstrap. Default 123.
fit: The fitted "nmfkc" object. Only needed when object does not carry X / Theta.
truncate: Logical; clip $Y^\ast$ at 0 (default TRUE).
epsilon, maxit: Convergence control of the re-fits, default 1e-8 and 20000 — tight, because a coefficient heading for the non-negativity boundary moves slowly and a loose tolerance biases anything that is compared with zero.
kkt.tol, kkt.ratio.min: Thresholds of the strict complementarity diagnostic (defaults 1e-3 and 5).
sep.tol: Row-normalised loading a column of $X$ must reach for that basis to count as anchored. Default 0.999. The zero threshold for the pure-column test on $\Theta$ is kkt.tol.
refit.args: Named list of extra arguments for the bootstrap re-fits. method and X.restriction are taken from the original fit automatically, so the replicates use the same estimator; anything else the original call set (X.init, penalties, ...) cannot be recovered from the fitted object and must be passed here. The function warns when it sees such arguments in the recorded call.
cores: Number of workers for the re-fits, default getOption("mc.cores", 1L) (sequential). The bootstrap multipliers are drawn before the loop and each re-fit is deterministic given its data, so the result does not depend on cores.

Value

An object of class "nmfkc.ar.latent.inference": a list with

coefficients: Data frame of the $DQ^2$ entries of the $G_d$, with columns Lag, Basis (effect at $t$), Covariate (cause at $t-d$), Estimate, BSE, CI_low, CI_high, p_value. Descriptive unless identified.approx.
G, G.se, G.ci.lower, G.ci.upper, G.p.value: The same quantities as lists of $Q\times Q$ matrices, one per lag.
spectral.radius: The radius of the $DQ\times DQ$ latent companion matrix — the same quantity nmfkc.ar.latent reports, not $\rho(\sum_d G_d)$, which differs when $D>1$. With spectral.radius.se, spectral.radius.ci.lower / .upper and spectral.radius.boot.mean.
prob.nonstationary: Bootstrap tail probability $P^\ast(\rho^\ast\ge1)$. Not a p-value for $H_0:\rho\ge1$ — see Details.
complex, complex.frac: Whether the point estimate has complex roots, and the fraction of replicates that do.
mu.y: With .ci.lower / .ci.upper; invariant, so inferential.
p.star: With .ci.lower / .ci.upper; descriptive unless identified.approx.
kkt: Strict-complementarity diagnostic: n.boundary, all.positive, ratio, and the two margins delta.dual and delta.prim.
identified.approx: Logical, an approximate diagnostic: TRUE for $Q=1$, for X.restriction = "fixed" (where $X$ is not estimated, so $T=I$), or when the per-column scale is fixed ("colSums" / "colSqSums") and the anchor-row and pure-column tests both pass at sep.tol / kkt.tol. It is not a certificate — the uniqueness argument needs exact anchors and exact zeros, and multiplicative updates leak a little. When FALSE the entry-level summaries and p.star are descriptive; $\rho$, the eigenvalues and $\bm\mu_y$ remain inferential either way.
identification: The two sides separately: anchor (per column of $X$), pure (per row of $\Theta$), ok, and the margins anchor.margin / pure.margin saying how far the fit is from an exact configuration (0 = exact).
refit.args: The arguments the bootstrap re-fits actually used, so the replicates can be checked against the original estimator.

Plus the bookkeeping entries wild.B, wild.B.requested, n.fail, truncate.rate, wild.unit, wild.dist, wild.level, dims.

Details

What is identified, and why that decides what can be inferred. The factorization admits $(X,\Theta)\to(XT,T^{-1}\Theta)$ for invertible $T$, so an entry of $\Theta$ means nothing across replicates; nmfkc.inference deals with this by conditioning on the estimated $X$. This function instead re-estimates the basis in every replicate, propagating its uncertainty. That splits the output in two:

Inferential. $\rho$, the eigenvalues and the observed long-run mean $\bm\mu_y=X\bm\mu_b$ are invariant under the whole family ($G_d\to T^{-1}G_dT$ leaves the spectrum alone), so their replicates need no alignment and their intervals mean what they say.
Conditional on identification. The entries $(G_d)_{qq'}$ and the composition $\bm p^\ast$ need $T$ to be a permutation, and three conditions together deliver that: the column scale fixed (X.restriction = "colSums"), an anchor row for every column of $X$, and a pure column for every row of $\Theta$. The anchor rows force $T\ge0$, the pure columns force $T^{-1}\ge0$, and a non-negative matrix with a non-negative inverse is monomial. Either one-sided condition alone is not enough: with $X=I_2$, $\Theta=\left(\begin{smallmatrix}2&3\\1&1\end{smallmatrix}\right)$ and $T=\left(\begin{smallmatrix}0.9&0.1\\0.1&0.9\end{smallmatrix}\right)$ the basis is perfectly separable and column-normalized, yet $XT\ge0$, $T^{-1}\Theta\ge0$, $X\Theta$ is unchanged, $T$ is no permutation, and $G$ differs between the two representatives — $\Theta$ has no pure column there. identified.approx reports the verdict and identification the two sides plus how far each is from exact; the function warns either way, naming the side that fails or the margins when it only passes approximately. The thresholds make this a diagnostic, not a certificate: the argument needs exact anchors and exact zeros, which multiplicative updates do not deliver. Replicates are permutation-aligned to the original basis regardless. X.restriction = "fixed" is the one clean case for $Q>1$ — the basis is not estimated, so $T=I$.

prob.nonstationary is not a unit-root test. It is the bootstrap tail probability $P^\ast(\rho^\ast\ge1)$, with replicates drawn around the estimated model. Nothing imposes $H_0:\rho=1$, so it does not have the interpretation of a p-value for that hypothesis; a proper test would have to resample under the null.

The zeros of $\Theta$ are only meaningful under misspecification. The strict-complementarity diagnostic reports the dual margin $\delta_{\mathrm{dual}}=\min|\Lambda^\ast|$ over the zero coefficients and the primal margin $\delta_{\mathrm{prim}}=\min\Theta$ over the positive ones. Both must be bounded away from zero for the support to be recovered exactly and for the active coefficients to attain the distribution that knowing the support would give. If the model is correctly specified, $\Xi_0=X\Theta$ makes $\Lambda^\ast=0$ exactly, so strict complementarity fails at every zero: $P(\hat\theta_k=0)$ then tends to a constant in $(0,1)$ (one half for an isolated zero) instead of to one, and no resampling scheme is consistent there. A small kkt$ratio is therefore not a numerical defect — it says the zeros are not estimating a sign restriction. What makes them estimate something is misspecification: they recover the coordinates whose unconstrained population coefficient is strictly negative.

Bootstrap validity at the boundary. A naive nonparametric bootstrap is inconsistent for a parameter on the boundary. The scheme used here resamples the residuals and re-imposes the non-negativity constraint in every re-fit, which is the practical remedy; a moving-block scheme would in addition preserve the serial dependence that the fixed design discards.

Bootstrap p-values invert the centred distribution. Every replicate of a non-negative coefficient is $\ge 0$, so comparing the raw replicates with zero would return $p = 0$ for every entry. The reported G.p.value is the basic-bootstrap two-sided $p$ obtained from $P^\ast(\hat G^\ast \ge 2\hat G)$.

Resampling scheme. With $\hat Y = X\Theta A$ and residuals $E=Y-\hat Y$, each replicate forms $Y^\ast=\hat Y+E\odot W$ and re-fits nmfkc at the same rank. $W$ carries one multiplier per time point (wild.unit = "column", the default, which preserves the contemporaneous covariance between series) or one per cell. Negative entries of $Y^\ast$ are truncated at 0 to keep the response non-negative; the truncation rate is reported, and a large one is a warning sign for the whole procedure.

This is a fixed-design bootstrap: the covariate matrix $A$ is held at the observed lags rather than rebuilt from $Y^\ast$, so the intervals are conditional on the observed past and do not propagate the randomness of the lagged design itself. They may therefore be somewhat optimistic.

References

Satoh, K. (2025). Applying non-negative matrix factorization with covariates to multivariate time series data as a vector autoregression model. Japanese Journal of Statistics and Data Science. arXiv:2501.17446. doi:10.1007/s42081-025-00314-0

Examples

# \donttest{
set.seed(1)
Y   <- matrix(abs(rnorm(4 * 40)) + 1, 4, 40)
ar  <- nmfkc.ar(Y, degree = 1)
fit <- nmfkc(ar$Y, ar$A, rank = 2, verbose = FALSE)
lat <- nmfkc.ar.latent(fit)
nmfkc.ar.latent.inference(lat, ar$Y, ar$A, wild.B = 50)
#> Warning: The factorization is not even approximately unique up to a permutation: basis 1, 2 has no anchor row in X (so T >= 0 is not forced); basis 1, 2 has no pure column in Theta (so T^-1 >= 0 is not forced). The entries of G_d and p.star are therefore DESCRIPTIVE; rho, the eigenvalues and mu.y are invariant under the whole family and remain the inferential targets.
#> Bootstrap inference for the latent VAR
#>   B = 50 (0 failed), unit = "column", dist = "rademacher", truncation rate = 0.32%
#> 
#> Latent transition matrices  G_d = Theta_d %*% X   (row = effect at t, column = cause at t-d)
#>   NOTE: the anchor-row / pure-column diagnostic FAILS, so these
#>        entries (and p.star) are DESCRIPTIVE. rho, the eigenvalues
#>        and mu.y are invariant and inferential.
#> 
#> lag1:
#>        Basis1 Basis2
#> Basis1  0.189  0.140
#> Basis2  0.085  0.083
#> SE:
#>        Basis1 Basis2
#> Basis1  0.152  0.141
#> Basis2  0.172  0.132
#> p-value:
#>        Basis1 Basis2
#> Basis1    0.2   0.28
#> Basis2    1.0   0.64
#> 
#> rho(G) = 0.2573   SE = 0.1484   95% CI = [0.0707, 0.5778]
#>   bootstrap mean = 0.3209  (bias +0.0636)
#>   bootstrap tail P*(rho* >= 1) = 0.000
#>   (a tail probability around the fitted model, NOT a unit-root test: the
#>    replicates are not drawn under an imposed H0: rho = 1)
#> 
#> complex eigenvalues: point estimate no, in 0.0% of replicates
#> 
#> long-run composition p*:
#>        estimate  lower  upper
#> Basis1   0.4749 0.2466 0.7462
#> Basis2   0.5251 0.2538 0.7534
# }