Bayesian Outlier Detection for Matrix-variate Models

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

86,624 characters · 13 sections · 53 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Bayesian Outlier Detection for Matrix--variate Models

abstractAnomalies in economic and financial data -- often linked to rare yet impactful events -- are of theoretical interest, but can also severely distort inference. Although outlier-robust methodologies can be used, many researchers prefer pre-processing strategies that remove outliers. In this work, an efficient sequential Bayesian framework is proposed for outlier detection based on the predictive Bayes Factor (BF). The proposed method is specifically designed for large, multidimensional datasets and extends univariate Bayesian model outlier detection procedures to the matrix-variate setting. Leveraging power-discounted priors, tractable predictive BF are obtained, thereby avoiding computationally intensive techniques. The BF finite sample distribution, the test critical region, and robust extensions of the test are introduced by exploiting the sampling variability. The framework supports online detection with analytical tractability, ensuring both accuracy and scalability. Its effectiveness is demonstrated through simulations, and three applications to reference datasets in macroeconomics and finance are provided.

Keywords: Bayesian Modelling, Bayes Factor, Matrix--variate, Sequential Model Assessment, Outliers.

Introduction

Robust methods are widely used to detect and address outliers -- observations that depart significantly from typical patterns due to measurement error, structural change, or data irregularities. In economics and finance, where such anomalies often signal impactful events like crises or policy shifts, proper identification is essential (atkinson1997detecting, GIORDANI2007112). Many researchers prefer to identify and remove outliers or windsorize the data (e.g., see malikov2020estimation, liao2022extrapolative, agarwal2024unobserved) rather than rely on robust methodologies as the latter are often more complex to implement and may not perform equally well across all scenarios (e.g., see zeng2021bayesian chang2024discussion). Instead, an extensive literature investigates the timing, magnitude, and probability of outliers following various robust model strategies such as change points (chopin2004bayesian, koop2007estimation, casini2024change), Markov switching (kole2023moments, casarin2024bayesian) and Bayesian nonparametrics (BASSETTI201449, BILLIO201997). In this paper, we propose a theoretically founded testing procedure for identifying frequency and periods of outlying observations. It is specifically designed for high-dimensional large datasets. Our test can be useful to the researcher for many purposes. When applied to the raw data, it can be used either to correctly choose a model class in the model specification stage or as a comparison with a model-based structural break identification. Additionally, it can be applied to the residuals of a given model to help identify sources of model mispecification.

Outlier detection is a statistical problem that has received considerable attention from the frequentist and Bayesian perspectives Bayarri2004. Despite their popularity, classical outlier detection procedures, such as Grubbs's test (Grubb50), are not designed for high-dimensional or structured data settings, such as multiple index panels or network-valued data. In such environments, observations often exhibit dependence, and univariate tests may yield misleading results.

Figure (ref) reports the outliers detected by Grubbs’ test across the three benchmark high-dimensional datasets we will consider in this work, i.e. an Inflation and Unemployment dataset (see Can09), an International Trade dataset (see Rose2004), and a Volatility Network dataset (see billio2018bayesianDT). The top panels display counts of individual outlying entries over time at different significance thresholds, together with the corresponding sample range (shaded areas). The bottom panels summarize the number of rows and columns containing at least one outlier. As illustrated in Figure (ref), applying standard outlier detection methods to the three high-dimensional datasets retrieves many outliers across different significance levels (top panels). This is to be expected, as several testing procedures that rely on independence assumptions tend to overlook interdependencies and therefore overflag outliers. When Bonferroni correction is applied, the number of detected outliers is substantially reduced due to the large number of series involved, highlighting the conservative nature of such inequality-based adjustment. Further results from Grubbs' test reveal structural dependencies and asymmetries between row- and column-wise detections (bottom panels), underscoring the need for joint testing. In addition, outlier detection for high-dimensional datasets calls for efficient and interpretable tools, which in turn motivate the development of tractable, iterative procedures for model fitting and prediction.

This paper addresses these issues and contributes to the Bayesian literature by extending sequential outlier detection procedures for univariate models to matrix-variate models. We follow a Bayesian approach as it naturally allows for sequential updating of the estimators involved in the outlier detection, and analytical results are derived to avoid computationally intensive approximations. The sequential nature of our test is well-suited for large and high-dimensional datasets as it allows reducing the computational cost, while also accommodating time variation in both the conditional mean and variance of the predictive distribution.

figure[figure omitted — 1,043 chars of source]

Within a Bayesian framework, a common approach is to assume that potential outliers arise from contaminating models distinct from the one generating the bulk of the data. Previous Bayesian works along these lines can be found in Box1968-gm,Guttman1973-bk, Abraham1979-ba, Guttman1978-be, pettit1985outliers,Pettit1992-bc,Pettit1990-zb,Verdinelli1991-cu,bayarri1994robust,Hoeting1996-ym. Nonetheless, these approaches are not well-suited for high-dimensional and large datasets since they are designed for a univariate setup.

We follow a Bayes Factor (BF) approach, which is also widely used in Bayesian analysis Schrider2016-nn,Ly2016-aw,Chen2018-ht,Li2021-pk, Stefan2019-nf, including selective inference Yekutieli2012-bd and outliers detection Bayarri2003-tb. The BF is also connected to Schwartz's criterion, also known as the Bayesian Information Criterion, which is another approach commonly used in the outlier detection literature. Due to the extreme sensitivity of the BF to the choice of the alternative distribution, robust methodologies have been developed Li2021-pk, Schad2022-oj. In this paper, we follow the power discounting approach of west1986bayesian and assume the alternative distribution is proportional to the principal distribution raised to a power.

The advantage of the perturbation approach based on power discounting is twofold. First, it requires estimating the model only under the null hypothesis, thus reducing computational cost. Second, the specification of the model under the alternative is not required, thus preserving the tractability of the BF testing procedure. We exploit these features to derive the finite-sample distribution of the predictive BF and to provide some frequentist validation of the testing procedure, without relying on asymptotic approximations. A critical region for the test is derived as an alternative to the traditional Jeffreys' scale of evidence, which has been criticized. See xian209 and referenced therein. We consider the Gaussian family since it is a standard assumption in many fields billio2018bayesianDT,guhaniyogi2017bayesian and extend the existing procedure for univariate Gaussian models west1986bayesian to the matrix--variate case, thus providing an original contribution to the expanding literature in this area landim2000dynamic,triantafyllopoulos2008missing,wang2009bayesian,CarvWest07DynMatNormGraph,Viroli11MatNorm,thompson2020classification,billio2018bayesianDT,tomarchio2022mixtures. The Gaussian assumption is not restrictive since we show that within a Bayesian framework the scale parameter can be integrated out of the likelihood, thus returning a Student-t predictive distribution, which can account for heavy tails.

This paper presents not only new results for detecting outliers in the matrix case but also uncovers novel results for the univariate Gaussian model. Furthermore, some pre-existing results for the univariate model are recovered as special cases. We also prove analytically that the outcome of the BF procedure for detecting outliers heavily depends on the choice of the discounting factor. Thus, we propose two robust BFs and a BF calibration procedure to alleviate the problem while maintaining a certain degree of tractability. Through some simulation experiments, we investigate the properties of a testing procedure under various outlier generation settings.

The impact of outliers on empirical economic analysis has gained importance in the aftermath of major global disruptions such as the 2008–2009 financial crisis and the COVID-19 pandemic. These episodes encouraged both researchers and official agencies to develop guidelines for outlier detection and adjusting for outliers in macroeconomic and financial data. In response to this growing need, we illustrate the effectiveness of our proposed sequential outlier detection procedure across the three aforementioned representative economic datasets: (i) a panel of inflation and unemployment indicators for European countries, (ii) a dynamic network of international trade flows, and (iii) a dynamic network of financial market volatilities.

The paper is organized as follows. Section (ref) introduces the outlier detection procedure based on BF and presents the main results. Section (ref) provides some analytical results on the properties of the test and a simulation study of the procedures. Section (ref) presents the three real-data illustrations. Section (ref) concludes.

Bayesian Outlier Detection

Bayes Factor

Consider a sequence of observations $\boldsymbol{Y}_t$, $t=1,2,\ldots$ with $\boldsymbol{Y}_{t}\in\mathcal{Y}$ where $\mathcal{Y}$ is a possibly multidimensional sample space and $t$ is a time index. In the following, boldfaced symbols represent vectors or matrices. We assume the information available at time $t$ is given by the collection of past observations $\mathcal{D}_t=\{\boldsymbol{Y}_1,\ldots,\boldsymbol{Y}_t\}$. Given the parameter and past observations, we assume the conditional sampling distribution does not depend on $\mathcal{D}_{t-1}$ and belongs to a parametric family with unknown parameter $\boldsymbol{\theta}\in \Theta$. The parameter space $\Theta$ is endowed with a prior distribution $p(\boldsymbol{\theta})$.

Our outlier detection procedure is applied sequentially over time to reduce computational cost in large datasets and capture time variations in the moments. We assume at time $t-1$ a posterior distribution for $\boldsymbol{\theta}$ is formed given $\mathcal{D}_{t-1}$ with density given by $p(\boldsymbol{\theta}|\mathcal{D}_{t-1})\propto p(\boldsymbol{\theta})p(\boldsymbol{Y}_1|\boldsymbol{\theta})\cdot\ldots\cdot p(\boldsymbol{Y}_{t-1}|\boldsymbol{\theta})$. We denote with $\boldsymbol{\theta}_t$ the random variable $\boldsymbol{\theta}|\mathcal{D}_{t-1}$ and assume its prior distribution is $p(\boldsymbol{\theta}_t |\mathcal{D}_{t-1})$. Assuming a model $g(\boldsymbol{Y}_t|\boldsymbol{\theta}_t)$ for $\boldsymbol{Y}_t$, the marginal predictive distribution for $\boldsymbol{Y}_t$ has density: $p(\boldsymbol{Y}_t |\mathcal{D}_{t-1})=\int_{\Theta}g(\boldsymbol{Y}_t|\boldsymbol{\theta}_t )p(\boldsymbol{\theta}_t |\mathcal{D}_{t-1}) d\boldsymbol{\theta}_t,\quad \boldsymbol{Y}_t\in\mathcal{Y}$. The distribution $p(\boldsymbol{Y}_t |\mathcal{D}_{t-1})$ naturally provides a measure of the predictive ability of the model $g(\boldsymbol{Y}_t|\boldsymbol{\theta}_t)$. In the following, we assume the predictive distribution belongs to the same distribution family as the sampling distribution, that is, $g(\boldsymbol{Y}_t|\boldsymbol{\theta}_t)=p(\boldsymbol{Y}_t|\boldsymbol{\theta}_t)$. In hypothesis testing, a model for the alternative hypothesis $\mathcal{H}_{1}$ is assumed, $p_A(\boldsymbol{Y}_t |\mathcal{D}_{t-1})$ and the prior predictive distribution $p_A(\boldsymbol{Y}_t |\mathcal{D}_{t-1})=\int_{\Theta}p(\boldsymbol{Y}_t |\boldsymbol{\theta}_t )p_A(\boldsymbol{\theta}_t |\mathcal{D}_{t-1}) d\boldsymbol{\theta}_t,\quad \boldsymbol{Y}_t\in\mathcal{Y}$ is obtained as the marginal distribution with respect to the prior distribution at time $t$ under the alternative hypothesis. The Bayes Factor is defined as the ratio of the two marginal distributions:

equation[equation omitted — 374 chars of source]

and yields an optimal decision in an inference problem under a 0-1 loss function. If $H_t>1$, the null hypothesis $\mathcal{H}_{0}$ is not rejected, then we conclude there is evidence from the observation $\boldsymbol{Y}_t$ in favor of the null hypothesis. This paper considers the BF as a testing procedure for model failure, which includes sensitivity to structural changes and outliers. To define the distribution under the alternative, we follow an approach based on a parametric perturbation $p_{A}(\boldsymbol{\theta}_t|\mathcal{D}_{t-1})= G(p(\boldsymbol{\theta}_t\vert \mathcal{D}_{t-1}),\alpha_t)$ of the posterior predictive distribution, where $\alpha_t$ is a possibly time-varying perturbation parameter and $G(\cdot,\cdot)$ is a perturbation function from $[0,1]\times[0,1]$ onto $[0,1]$. The main advantages of the perturbation approach are two. First, only the estimation of the model under the null is required, which substantially reduces the computational cost in large datasets and high-dimensional models. Second, the specification of the model under the alternative is not required, which usually demands rather involved models and a high computational cost. On the contrary, the perturbation approach preserves the tractability of the model under the null and provides an effective strategy for capturing deviations from the null hypothesis.

Examples of perturbation functions can be derived from the approaches based on mixtures of distributions, $p_{A}(\boldsymbol{\theta}_t\vert \mathcal{D}_{t-1})=\alpha_t p(\boldsymbol{\theta}_t)+(1-\alpha_t)p(\boldsymbol{\theta}_t\vert \mathcal{D}_{t-1})$ robert202250, and on distortion of probability measures, $p_{A}(\boldsymbol{\theta}_t\vert \mathcal{D}_{t-1})= g(\int_{-\infty}^{\boldsymbol{\theta}_t}p(u\vert \mathcal{D}_{t-1})du,\alpha_t)p(\boldsymbol{\theta}_t\vert \mathcal{D}_{t-1})$ where $g$ is the partial derivative of $G$ with respect to its first argument.

In this paper, we assume the following perturbation function $G(p,\alpha_t)=p^{\alpha_t} C(\alpha_t)$ where $C(\alpha_t)=\left(\int_{\Theta}p(\boldsymbol{\theta}_t |\mathcal{D}_{t-1})^{\alpha_t}d\boldsymbol{\theta}_t\right)^{-1}$ is the inverse normalizing constant of the density under the alternative hypothesis. The constant satisfies $C(\alpha_t)\rightarrow 1$ as $\alpha_t\rightarrow 1^{-}$ and $C(\alpha_t)\rightarrow\lambda(\Theta)^{-1}$ as $\alpha_t\rightarrow 0^{+}$, where $\lambda(\Theta)<\infty$ denotes the Lebesgue measure of $\Theta$. If $\lambda(\Theta)$ is unbounded, then to prevent $C(\alpha_t)\rightarrow 0$ as $\alpha_t\rightarrow 0^{+}$, some restrictions can be introduced such as $\alpha_t\in(\underline{\alpha},\bar{\alpha})\subset(0,1)$, $\underline{\alpha}>0$. These limits are relevant in the analysis of the BF, denoted as $H_{t}(\alpha_t)$ in what follows. This choice of the calibration function returns the power discounting approach proposed in the seminal paper west1986bayesian for eliciting the prior under the alternative hypothesis

equation[equation omitted — 142 chars of source]

where the discounting parameter $\alpha_t$ takes values in the interval $(0,1)$. An advantage of the power discounting approach is that it returns the popular variance contamination model for outlier detection Pettit1992-bc,weiss1997bayesian,Bayarri2003-tb,page2011bayesian,tomarchio2022mixtures and does not require the estimation of the contamination parameter, which would be difficult, since very little information is available about the outlier generating process. See also van2021bayesian for a comparison between Bayesian contamination procedures and alternative approaches.

We assume a Gaussian likelihood and a conjugate prior distribution to preserve some analytical tractability and extend the univariate discounting approach to a matrix--variate setting. We will prove that the BF is generally well-defined and bounded by a function of $\alpha_t$. However, the results of the hypothesis testing can be significantly influenced by the value of $\alpha_t$, and the BF can become unbounded in the limit for $\alpha_t\rightarrow 0^{+}$. Consequently, the discounting parameter should be chosen carefully. The new results for the matrix setting readily apply to the previous univariate and multivariate approaches.

Robust Bayes Factors

As one can expect, the outcome of a decision based on the BF may depend on the value of $\alpha_t$. See also the discussion in west1986bayesian. We prove that for a given level $H^{\mathbin{\vcenter{\hbox{\scalebox{.5}{$\bullet$}}}}}$ of the threshold (which is assumed equal to one in the following), for some values of $\alpha_t$, the evidence is against the null hypothesis, and for some others, it is against the alternative. In the outlier detection setting, this means that a new observation $Y_t$ may be considered an outlier or not based on the value of $\alpha_t$. The following result illustrates the indeterminacy of the outcome of the BF procedure.

propositionAssume the likelihood satisfies \begin{equation} \int_{\Theta}p(\boldsymbol{Y}_t|\boldsymbol{\theta}_t)(\lambda(\Theta)p(\boldsymbol{\theta}_t|\mathcal{D}_{t-1})-1)\lambda(d\boldsymbol{\theta}_t)>0 \end{equation} with respect to the reference measure $\lambda$, and there exists a unique $\alpha_{t}^{\mathbin{\vcenter{\hbox{\scalebox{.5}{$\bullet$}}}}}$ such that $q_1(\alpha_{t})-q_2(\alpha_{t})<0$ for $\alpha_t<\alpha_{t}^{\mathbin{\vcenter{\hbox{\scalebox{.5}{$\bullet$}}}}}$ and $q_1(\alpha_{t})-q_2(\alpha_{t})>0$ for $\alpha_t>\alpha_{t}^{\mathbin{\vcenter{\hbox{\scalebox{.5}{$\bullet$}}}}}$ and $q_1(\alpha_{t}^{\mathbin{\vcenter{\hbox{\scalebox{.5}{$\bullet$}}}}})-q_2(\alpha_{t}^{\mathbin{\vcenter{\hbox{\scalebox{.5}{$\bullet$}}}}})=0$, where \begin{eqnarray} &&q_1(\alpha_t)=\int_{\Theta}p(\boldsymbol{Y}_t|\boldsymbol{\theta}_t)p(\boldsymbol{\theta}_t|\mathcal{D}_{t-1})^{\alpha_t}\lambda(d\boldsymbol{\theta}_t) \int_{\Theta}p(\boldsymbol{\theta}_t|\mathcal{D}_{t-1})^{\alpha_t}\log p(\boldsymbol{\theta}_t|\mathcal{D}_{t-1})\lambda(d\boldsymbol{\theta}_t)\\ &&q_2(\alpha_t)=\int_{\Theta}p(\boldsymbol{\theta}_t|\mathcal{D}_{t-1})^{\alpha_t}\lambda(d\boldsymbol{\theta}_t)\int_{\Theta}p(\boldsymbol{Y}_t|\boldsymbol{\theta}_t)p(\boldsymbol{\theta}_t|\mathcal{D}_{t-1})^{\alpha_t}\log p(\boldsymbol{\theta}_t|\mathcal{D}_{t-1})\lambda(d\boldsymbol{\theta}_t). \end{eqnarray} Then there exists a partition $(U_1,U_2)$ of the unit interval such that $H_t(\alpha_t)>1$ for $\alpha_{t}\in U_1$ and $H_t(\alpha_t)<1$ for $\alpha_{t}\in U_2$.

Proposition (ref) establishes the conditions such that different values of $ \alpha_t $ may cause the BF $ H_t(\alpha_t) $ to uniquely cross the threshold value $ H_t(\alpha^{\mathbin{\vcenter{\hbox{\scalebox{.5}{$\bullet$}}}}}_{t}) = 1 $ at $ \alpha_{t}^{\mathbin{\vcenter{\hbox{\scalebox{.5}{$\bullet$}}}}} $. It can be seen from the proof that for $\lambda(\Theta)$ unbounded, the evidence against the presence of an outlier is negative for $\alpha_t\rightarrow0^{+}$, consistently with the Jeffreys--Lindley Paradox Robert_2014, bernardo2009bayesian. See also robert2009harold and wag2023 for a review with a historical perspective. In addition, we shall notice that, in the Bayesian current practice, the BF is usually interpreted following Jeffrey's scale of evidence against the null hypothesis kass1995bayes, which suggests the evidence is negative for $1/H_t<1$, not worth more than a bare mention for $1<1/H_t<10^{1/2}$, substantial for $10^{1/2}<1/H_t<10$, strong for $10<1/H_t<10^{3/2}$, very strong for $10^{3/2}<1/H_t<10^2$ and decisive for $1/H_t>10^2$. The following example illustrates the result given in the previous proposition and shows that an observation can be regarded as an outlier depending on the value of $\alpha_t$. The illustration is general since a stepwise uniform prior is assumed, and the differentiability of the prior is not required to find the threshold value of $\alpha_t$.

rem[Univariate Gaussian model] Assume a Gaussian likelihood with location $\theta_t$ and variance $\sigma^2$, that is $Y_t \sim\mathcal{N}(\theta_t,\sigma^2)$ independent for $t=1,\ldots,T$. Let us assume $\sigma^2$ is known and $\theta_t$ follows a prior distribution $p(\theta_t\vert \mathcal{D}_{t-1})$ which is absolutely continuous with respect to the Lebesgue measure $\lambda$, and has density function \begin{equation} p(\theta_t\vert \mathcal{D}_{t-1})=\sum_{j=1}^{n}g_{t,j}\mathbb{I}(\theta_t\in \Theta_j), \end{equation} where $\{\Theta_j,j=1,\ldots,n\}$ and $\Theta=\Theta_1\cup\ldots\cup\Theta_n$ with $n>1$ and such that $g_{t,1}\lambda(\Theta_1)+\ldots+g_{t,n}\lambda(\Theta_n)=1$. In Table (ref), we report three cases in which an observation $Y_t$ can be regarded as an outlier depending on the value of $\alpha_t$. The table reports together with $Y_t$, the value $\alpha^{\mathbin{\vcenter{\hbox{\scalebox{.5}{$\bullet$}}}}}_t$ such that $H_t(\alpha^{\mathbin{\vcenter{\hbox{\scalebox{.5}{$\bullet$}}}}}_t) = 1$, the step-function prior specification with possibly disconnected support given by the union of $\Theta_j=(\underline{\theta}_{t,j},\overline{\theta}_{t,j})$, and density $g_{t,j}$ for $j \in 1, \ldots, 3$. The table reports also the likelihood values $p_{t,j}=p(Y_t|\theta_{t,j})$ with $\theta_{t,j} \in \Theta_j$. \begin{table}[h!] {2pt} \begin{tabular}{ccc|cccccc|ccc|ccc} \hline $\alpha_{t}^{\mathbin{\vcenter{\hbox{\scalebox{.5}{$\bullet$}}}}}$& $Y_t$ &$\sigma$&$\underline{\theta}_{t,1}$ & $\overline{\theta}_{t,1}$ & $\underline{\theta}_{t,2}$ & $\overline{\theta}_{t,2}$ &$\underline{\theta}_{t,3}$ &$\overline{\theta}_{t,3}$&$g_{t,1}$&$g_{t,2}$&$g_{t,3}$&$p_{t,1}$&$p_{t,2}$&$p_{t,3}$\\ \hline 0.054&-4.422&1.120&0.471&0.565&1.487&1.507&2.863&5.863&0.440&50.000&0.001&0.180 & 0.001&0.001\\ 0.076&-5.690&1.054&-1.112&1.689&2.422&2.588&4.780&4.800&0.002&0.492&45.645&0.002 & 0.153&0.002\\ 0.412&-8.723&1.939&0.996&1.409&2.628&4.209&5.498&5.549&0.934&0.023&11.510&0.146 & 0.019&0.041\\ \hline \end{tabular} \caption{Prior hyper-parameters $\underline{\theta}_{t,j}$, $\overline{\theta}_{t,j}$ and $g_{t,j}$ and the likelihood values $p_{t,j}=p(Y_t|\theta_{t,j})$, with $j=1,\ldots,3$.} \end{table} The marginal likelihood under the alternative hypothesis is derived as follows. By the mean-value theorem for integrals on a bounded domain, there exists $\theta_{t,j}\in\Theta_j$ such that \begin{equation} \int_{\Theta}p(Y_t|\theta_t)p_A(\theta_t|\mathcal{D}_{t-1})d\theta_t=C(\alpha_t)\sum_{j=1}^{n}\int_{\Theta_j} p(Y_t\vert \theta_t)g_{t,j}^{\alpha_t}d\theta_t =C(\alpha_t)\sum_{j=1}^{n}g_{t,j}^{\alpha_t}p_{t,j}\lambda(\Theta_j), \end{equation} where $p_{t,j}=p(Y_t|\theta_{t,j})$ and $C(\alpha_t)=g_{t,1}^{\alpha_t}+g_{t,2}^{\alpha_t}+g_{t,3}^{\alpha_t}$ is the normalizing constant. This expression can is used to show that the conditions in Prop. (ref) are satisfied for the cases provided in Table (ref). The values of $p_{t,j}$ show that the prior predictive assigns a non--negligible probability to the observations considered in the example. Figure (ref) shows that the observation $Y_t$ is an outlier for values of $\alpha_t < \alpha_t^{\mathbin{\vcenter{\hbox{\scalebox{.5}{$\bullet$}}}}}$ ($H_t(\alpha_t)>1$), while $Y_t$ is considered an outlier for $\alpha_t>\alpha_t^{\mathbin{\vcenter{\hbox{\scalebox{.5}{$\bullet$}}}}}$ ($H_t(\alpha_t)<1$). Also, for some observations and settings, comparing the BF with the first threshold of Jeffrey's scale, $10^{1/2}$, returns substantial evidence against the null (first and second panel) depending on $\alpha_t$. The third example shows that the farther the observed value from $\Theta$, e.g. $Y_t=-8.723$, the larger $\alpha_t$ needs to be in order to make the presence of an outlier a convincing hypothesis, i.e. $\alpha_{t}^{\mathbin{\vcenter{\hbox{\scalebox{.5}{$\bullet$}}}}}=0.412$.

figure[figure omitted — 696 chars of source]

Since the outcome of the testing procedure strongly depends on the choice of the discounting parameter $\alpha_t$ and since different thresholds for $H_t$ can be chosen, we propose three alternative Bayesian decision rules to conduct outlier detection. These rules are robust and exploit the variability of the BF due to the discounting parameter and the sampling distribution.

A first decision rule is the Minimum BF (MBF), i.e. the smallest possible BF within the class of alternative distributions specified in Eq. (ref), that is $$H_t^{(MBF)}=\inf_{\alpha_t\in(\underline{\alpha},\overline{\alpha})\subset (0,1)} H_t(\alpha_t).$$ The MBF has been used in other testing problems such as unit root testing berger1994noninformative, correlation testing Chen03062021 and reverse-Bayes procedures pawel2022sceptical. See also HeldOtt2017 for a review of minimum BFs and a comparison with standard BF. The rationale behind this decision rule is that the evidence for the null is at least the MBF. In Remark (ref), the three settings return an MBF below one. However, the MBF is not too far below the threshold in the first setting compared to the other two scenarios. These considerations motivate the need for alternative decision rules.

We consider the Integrated BF (IBF) as a second decision rule. With this rule, one assumes a prior distribution for the discounting parameter $\alpha_t$ and averages the BF over all possible discounting values, that is $$H_t^{(IBF)}=\int_{\underline{\alpha}}^{\bar{\alpha}} H_t(\alpha_t)\pi_t(\alpha_t)d\alpha_t-1,$$ where $\pi_t(\alpha)$ is a suitable probability density function for $\alpha_t$ with support $(\underline{\alpha},\bar{\alpha})\subset(0,1)$. The choice of $\pi_t$ is crucial to achieve a well-defined IBF in the matrix variate case.

Since the BF can be bounded from above under mild regularity conditions, the IBF can be modified to account for the relative magnitude of the evidence in favor of the null. We thus introduce the Normalized IBF (NIBF) defined as $$H_t^{(NIBF)}=\left(\int_{\underline{\alpha}}^{\bar{\alpha}} H_t(\alpha_t)\pi_t(\alpha_t)d\alpha_t-1\right)\left(\int_{\underline{\alpha}}^{\bar{\alpha}} \kappa_t(\alpha_t)\pi_t(\alpha_t)d\alpha_t-1\right)^{-1},$$ where $\kappa_t(\alpha_t) < \infty$ is an upper bound for $H_t(\alpha_t)$. The upper bound will be used later in this paper to show the BF integrability with respect to the discounting parameter.

IBF and NIBF do not account for the sampling variability in hypothesis testing. For this reason, we propose a predictive BF approach to incorporate such variability. As $Y_t$ is not observed at time $t-1$, the predictive BF is a random variable whose marginal distribution can be used to derive a calibrated value for the discounting parameter $\alpha_t$. With our approach, we account for the sampling variability of the predictive BF and find the analytical distributions $F_{j,t}$ of the random BF, $H_{t}(\alpha_t)$, under the null hypothesis $\mathcal{H}_0$ of the absence of outliers ($j=0$) and the alternative hypothesis $\mathcal{H}_1$ of the presence of outliers ($j=1$). The calibrated predictive BF is derived through the following steps.

enumerate[1.] • Derive the threshold function $\underline{h}(\alpha_t)$, $\alpha_{t}\in(\underline{\alpha},\bar{\alpha})$ for rejecting the null $\mathcal{H}_0$, such that $\mathbb{P}(\{H_t<\underline{h}(\alpha_t)\})=\tau$ under the null for a given test size $\tau$. The probability is evaluated using the $F_{0,t}$ distribution. • Determine the discounting parameter $\alpha_t^*\in(\underline{\alpha},\bar{\alpha})$ such that for a given power level $\beta$ $\mathbb{P}(\{H_t>\bar{h}(\alpha_t^*)\})=1-\beta$ under the alternative, that is: $\beta = 1 - F_{1,t}(\bar{h}(\alpha_t^*))$, $\bar{h}(\alpha_t^*) = 2 - \underline{h}(\alpha_t^*)$, where the constraint comes from the threshold symmetric assumption. • Define the inconclusive interval $C_H = (\underline{h}(\alpha_t^*), \bar{h}(\alpha_t^*))$ where the BF is statistically equal to one and define the calibrated BF as $H_t^{*}=H_t(\alpha_t^{*})$.

The conclusion of the test procedure is to reject the null hypothesis $\mathcal{H}_0$ when $H_t<\underline{h}$, accepting it when $H_t>\bar{h}$ and randomizing when $H_t\in C_{H}$. As randomization is usually not appealing in some applications, an alternative procedure can be used where only one threshold $\bar{h}(\alpha_t^*)=\underline{h}(\alpha_t^*)<1$ is considered to define a critical region $(0,\underline{h}(\alpha_t^*))$. Our procedure for the calibrated value of $\alpha_t$ is similar in spirit to the ones proposed in weiss1997bayesian and pawel2025closed for determining the optimal sample size. In addition, another procedure, which accounts for the distribution of $\alpha_t$, can be defined by minimizing jointly the first and second type error probabilities.

Matrix Normal Outlier Detection

A Bayesian Matrix Normal Model

The BF and its properties are derived under the following assumptions. See the Appendix (ref) for some background on the matrix distributions and proofs. The normal assumption is standard in outlier detection as it guarantees some finite-sample analytical results.

assumption$\mathbf{Y}=\{\boldsymbol{Y}_{j}\left(p\times n\right),\ j = 1,\ldots,t - 1\}$, are i.i.d. matrix normal variables with distribution $\mathcal{MN}_{p,n}(\boldsymbol{B},\boldsymbol{\Sigma}_L,\boldsymbol{V})$ with $\boldsymbol{B}\in\mathbb{R}^{p\times n}$ $\boldsymbol{\Sigma}_{L}\in\mathbb{R}^{p\times p}_{+} $ and $\boldsymbol{V}\in\mathbb{R}^{n\times n}_{+}$.

The information set at time $T=t-1$ is given by the sigma-algebra $\mathcal{D}_{t-1}$ generated by the elements of $\mathbf{Y}$. To achieve analytical tractability, and similarly to the univariate setting of west1986bayesian, we assume that $\boldsymbol{\Sigma}_{L}$ is $\mathcal{D}_{t-1}$-measurable and $\boldsymbol{B}$ has a conjugate prior. We shall emphasize that the Gaussian assumption is not restrictive, as within a Bayesian framework, the scale parameters can be integrated out of the likelihood, resulting in a Student-t distribution that can account for heavy tails.

assumption$\boldsymbol{B}$ has a conjugate matrix normal prior, i.e. $\boldsymbol{B}\sim \mathcal{N}_{p,n}\left( \boldsymbol{M},\boldsymbol{\Sigma}_{P},\boldsymbol{V} \right)$ and $\boldsymbol{\Sigma}_{L}\otimes \boldsymbol{V}$ is given.

The assumption $\boldsymbol{\Sigma}_{P} = \boldsymbol{\Sigma}_{L}/\varphi$, with $\varphi > 0$ is common in Gaussian models zellner1986assessing. Some analytical results can also be obtained when the covariance $\boldsymbol{V}$ is estimated following a Bayesian procedure. When $\boldsymbol{V}$ is unknown, a conjugate Normal-Inverse Wishart prior for $\boldsymbol{B}$ and $\boldsymbol{V}$ is assumed.

assumption$\boldsymbol{B}$ and $\boldsymbol{\Sigma}_L$ have a conjugate inverse Wishart matrix normal prior, i.e. $\boldsymbol{B}|\boldsymbol{V}\sim \mathcal{N}_{p,n}\left( \boldsymbol{M},\boldsymbol{\Sigma}_{L}/\varphi,\boldsymbol{V}/\rho \right)$ and $\boldsymbol{V}\sim\mathcal{I}\mathcal{W}_n(\boldsymbol{\Psi},m)$ and $\boldsymbol{\Sigma}_{L}$ is given.

Following the notation in the previous section, $\boldsymbol{\theta}_{t}$ corresponds to $\boldsymbol{B}_t=\boldsymbol{B}|\mathcal{D}_{t-1}$ and $(\boldsymbol{B}_t,\boldsymbol{V}_t)=(\boldsymbol{B},\boldsymbol{V})|\mathcal{D}_{t-1}$ for the $\boldsymbol{V}$ known and $\boldsymbol{V}$ unknown cases, respectively.

Proposition (ref) provides the analytical expression of the Bayes Factor derived under the three assumptions outlined above.

proposition[BF] \begin{enumerate}[i)] • Under the Assumptions (ref) and (ref), $\boldsymbol{Y}_t|\mathcal{D}_{t - 1}\sim \mathcal{N}_{p,n}\left(\boldsymbol{M}_{*},\boldsymbol{\Sigma}_{d}, \boldsymbol{V} \right)$ under the null hypothesis, and $\boldsymbol{Y}_t|\mathcal{D}_{t - 1}\sim \mathcal{N}_{p,n}\left(\boldsymbol{M}_{*},\boldsymbol{\Sigma}_{A,d}, \boldsymbol{V} \right)$ under the alternative and the Bayes Factor $H_t$ is given by \begin{equation} H_{t}(\alpha_{t}) = \frac{\left| \boldsymbol{\Sigma}_{A,d} \right|^{\frac{n}{2}}}{\left| \boldsymbol{\Sigma}_{d} \right|^{\frac{n}{2}}}\exp\left\{ - \frac{1}{2}tr\left\lbrack \left( \boldsymbol{\Sigma}_{d}^{- 1} - \boldsymbol{\Sigma}_{A,d}^{- 1} \right)\left( \boldsymbol{Y}_t - \boldsymbol{M}_{*} \right)\boldsymbol{V}^{- 1}\left( \boldsymbol{Y}_t - \boldsymbol{M}_{*} \right)^{'} \right\rbrack \right\}, \end{equation} where ${\boldsymbol{\Sigma}_{d} = \boldsymbol{\Sigma}}_{L} + \boldsymbol{\Sigma}_{*}$, $\boldsymbol{\Sigma}_{A,d} = \boldsymbol{\Sigma}_{L} + \boldsymbol{\Sigma}_{A,*}$, $\boldsymbol{\Sigma}_{A,*} = \boldsymbol{\Sigma}_{*}/\alpha_{t}$ and the posterior parameters $\boldsymbol{M}_{\ast}$ and $\boldsymbol{\Sigma}_{*}$ are given in Eq. (ref) and (ref). • Under the the Assumptions (ref) and (ref), $\boldsymbol{Y}_t|\mathcal{D}_{t - 1}\sim\mathcal{T}_{p,n}\left(m_{*}-2n, \boldsymbol{M}_{*}, \boldsymbol{\Sigma}_L , \boldsymbol{L}_{*} \right)$ under the null and $\boldsymbol{Y}_t|\mathcal{D}_{t - 1} \sim \mathcal{T}_{p,n}\left(m_{A, *}-2n, \boldsymbol{M}_{*}, \boldsymbol{\Sigma}_L , \boldsymbol{L}_{A, *}\right)$ under the alternative where $\mathcal{T}_{p,n}( \nu, \boldsymbol{M},\boldsymbol{\Sigma}, \boldsymbol{\Omega})$ denotes the matrix Student-t distribution, $\boldsymbol{L}_{*} = \boldsymbol{\Psi}_{*}k_d/k_{*}$, and $\boldsymbol{L}_{A,*} = \boldsymbol{\Psi}_{A,*}k_d/k_{A,*}$. The Bayes Factor $H_t$ is given by \begin{equation} H_t(\alpha_t)=G\frac{\left|\boldsymbol{\Psi}_{A,*}+\frac{\alpha_t k_{*}}{\left(\alpha_t k_{*}+1\right)}\left(\boldsymbol{M}_{*}-\boldsymbol{Y}_t\right)^{'}\boldsymbol{\Sigma}_L^{-1}\left(\boldsymbol{M}_{*}-\boldsymbol{Y}_t\right)\right|^{\frac{m_{A,d}-n-1}{2}}}{\left|\boldsymbol{\Psi}_{*}+k_{*}\frac{\left(\boldsymbol{M}_{*}-\boldsymbol{Y}_t\right)^{'}\boldsymbol{\Sigma}^{-1}\left(\boldsymbol{M}_{*}-\boldsymbol{Y}_t\right)}{\left(k_{*}+1\right)}\right|^{\frac{m_{d}-n-1}{2}}}, \end{equation} where \begin{equation} G=\frac{\left|\boldsymbol{\Psi}_{*}\right|^{\frac{m_{*}-n-1}{2}}k_{*}^\frac{np}{2}\Gamma_{n}\left(\frac{m_{d}-n-1}{2}\right)k_{A,d}^\frac{np}{2}\Gamma_{n}\left(\frac{m_{A,*}-n-1}{2}\right)}{k_{d}^\frac{np}{2}\Gamma_{n}\left(\frac{m_{*}-n-1}{2}\right)\left|\boldsymbol{\Psi}_{A,*}\right|^{\frac{m_{A,*}-n-1}{2}}k_{A,*}^\frac{np}{2}\Gamma_{n}\left(\frac{m_{A,d}-n-1}{2}\right)}\nonumber\\ \end{equation} and $k_{d}=k_{*}+1$, $k_{A,*}=\alpha_t k_{*}$, $ k_{A,d}=\alpha_t k_{*}+1$, $m_{d}=m_{*}+p$, $m_{A,*}=\alpha_t( m_{*}+p)-p$, $m_{A,d}=m_{A,*}+p$, and $m_{*}, \boldsymbol{M}_{*},\boldsymbol{\Psi}_{*}$ as in Eq. (ref). \end{enumerate}

rem[Interpretability] In the case $\boldsymbol{V}$ is known, it can be easily shown that the discounting approach provides a tractable prior distribution under the alternative. The approach is equivalent to assuming, under the alternative hypothesis, the following hierarchical prior distribution $\boldsymbol{B}|\tilde{\boldsymbol{M}},\mathbf{Y}\sim\mathcal{N}_{p,n}(\tilde{\boldsymbol{M}},\boldsymbol{\Sigma}_{*},\boldsymbol{V})$ with $\tilde{\boldsymbol{M}}|\mathbf{Y}\sim\mathcal{N}_{p,n}(\boldsymbol{M}_{*},\boldsymbol{\Sigma}_{*}(1/\alpha_t-1),\boldsymbol{V})$, which is well-defined for $\alpha_t \in (0,1)$. This can be interpreted as a random perturbation of the posterior distribution, which inflates the variance of the location parameter distribution.

The normalizing constant of the alternative distribution and its properties can be derived from Prop. (ref) and are given in the following.

corollary\begin{enumerate} • Under Assumptions (ref) and (ref), the inverse normalizing constant is \begin{equation} C\left( \alpha_{t} \right)={\alpha_{t}^{\frac{np}{2}}(2\pi)}^{\frac{(\alpha_{t}-1)np}{2}}\left| \boldsymbol{\Sigma}_{*} \right|^{\frac{(\alpha_{t}-1)n}{2}}|\boldsymbol{V}|^{\frac{(\alpha_{t}-1)p}{2}} \end{equation} and satisfies $C(\alpha_t)\rightarrow 0$ as $\alpha_{t} \rightarrow 0^{+}$ and $C(\alpha_t)\rightarrow 1$ as $\alpha_t\rightarrow 1^{-}$. • Under Assumptions (ref) and (ref), the inverse normalizing constant is \begin{equation} C(\alpha_t) = \frac{{\left|\boldsymbol{\Psi}_{A,*}\right|^{\frac{m_{A,*}-n-1}{2}}}(2\pi)^{\frac{\alpha_{t}np}{2}}|\boldsymbol{\Sigma}_{*}|^{\frac{\alpha_{t}n}{2}}2^{\frac{\alpha_{t}(m_{*}-n-1)n}{2}}\Gamma_{n}\left(\frac{m_{*}-n-1}{2}\right)^{\alpha_{t}}}{(2\pi)^{\frac{np}{2}}|\boldsymbol{\Sigma}_{A,*}|^{\frac{n}{2}}2^{\frac{(m_{A,*}-n-1)n}{2}}\Gamma_{n}\left(\frac{m_{A,*}-n-1}{2}\right)\left|\boldsymbol{\Psi}_{*}\right|^\frac{\alpha_t (m_{*}-n-1)}{2}} \end{equation} \end{enumerate} and satisfies $C(\alpha_t)\rightarrow 0$ as $\alpha\rightarrow \underline{\alpha_t}^{+}$ and $C(\alpha_t)\rightarrow 1$ as $\alpha_t\rightarrow 1^{-}$, where $\underline{\alpha}=(p+2n)/m_d$.

rem[Univariate Gaussian model] Let $Y_{j}\sim\mathcal{N}(\theta,\sigma^2)$ i.i.d. for $j=1,\ldots,t-1$ and assume a conjugate prior for $\theta$, that is $\theta\sim\mathcal{N}(m,\sigma^2/\varphi)$ with $\varphi>0$. It is easy to show that the predictive under the alternative is given by $p_{A}(Y_{t}|\mathcal{D}_{t-1})=(2\pi\sigma^2_{\ast}/\alpha_t)^{1/2}\exp\{-(\theta-m_{\ast})^2\alpha_t/2\sigma^{2}_{\ast}\}$ with $m_{\ast}=\varphi m/(\varphi+t-1)+(t-1)/(\varphi+t-1)\bar{Y}$ with $\bar{Y} = (t-1)^{-1}\sum_{j = 1}^{t-1}Y_j$, $\sigma^2_{*}=\sigma^2/(\varphi+t-1)$ and the normalizing constant is $C(\alpha_t)=(2\pi\sigma^2_{*})^{(\alpha_t-1)/2}\alpha_t^{1/2}$. It follows that the BF is $H_t(\alpha_t)=\kappa_t(\alpha_t)\exp\{(Y_t-m_{\ast})^2 A_t^{-1}\}$, where \begin{equation} A_t=\frac{2\sigma^2_{\ast}(\varphi+t)(\alpha_t(\varphi+t-1)+1)}{(\alpha_t-1)(\varphi+t-1)}\quad\hbox{and}\quad \kappa_t(\alpha_t)=\left(\frac{\alpha_t(\varphi+t-1)+1}{\alpha_t(\varphi+t)}\right)^{1/2} \end{equation} naturally represents an upper bound of $H_{t}(\alpha_t)$. Intuitively, for small values of $\alpha_t$ the null hypothesis is not rejected, that is $H_t(\alpha_t)\rightarrow +\infty$ as $\alpha_t\rightarrow 0^+$, whereas for large values, the BF gets closer to one, that is $H_t(\alpha_t)\rightarrow 1$ as $\alpha_t\rightarrow 1^{-}$. Its derivative $H^{\prime}_{t}(\alpha_t)=0$ for $\alpha_t$ is equal to $\alpha_{t,0}=((Y_t-m_{\ast})^2/\sigma^2_{\ast}-1)^{-1}(\varphi+t-1)^{-1}$ provided $(Y_{t}-m_{\ast})^2/\sigma^2-1>(\varphi+t-1)^{-1}$.

Bayes Factor Properties

proposition[BF properties, known $\boldsymbol{V}$] Under Assumptions (ref) and (ref), the BF $H_{t}(\alpha_t)$ satisfies the following properties. \begin{enumerate}[i)] • There exists a function $\kappa_t\left( \alpha_{t} \right)= \left| \boldsymbol{\Sigma}_{A,d} \right|^{\frac{n}{2}}\left| \boldsymbol{\Sigma}_{d} \right|^{-\frac{n}{2}}\geq1$ independent of $\boldsymbol{Y}_t$ such that $H_t(\alpha_t)\leq \kappa_t(\alpha_t)$, $\kappa_t \rightarrow 1$ as $\alpha_{t} \rightarrow 1^{-}$, $\kappa_t \rightarrow \infty$ as $\alpha_{t} \rightarrow 0^{+}$ and $\kappa_t$ is decreasing in $\alpha_t$; • $H_{t}\left( \alpha_{t} \right) \rightarrow 1$ as $\alpha_{t} \rightarrow 1^{-}$ and $H_{t}\left( \alpha_{t} \right)\rightarrow\infty$ as $\alpha_{t} \rightarrow 0^{+}$; • $\partial_{\alpha_t} H_t(\alpha_t)= H_{t}\left( \alpha_{t} \right)(2\alpha_{t})^{-1}\text{tr}\left(-n\Tilde{\boldsymbol{B}}^{- 1}\boldsymbol{\Sigma}_{*} +\alpha_{t} \boldsymbol{\Upsilon}(\alpha_t) \Tilde{\boldsymbol{A}} + \alpha_{t}\left( {1 - \alpha}_{t} \right)\boldsymbol{\Upsilon}(\alpha_t)\boldsymbol{\Sigma}_{L}\Tilde{\boldsymbol{B}}^{-1}\Tilde{\boldsymbol{A}} \right)$ where we defined $\Tilde{\boldsymbol{A}} = \left( \boldsymbol{Y}_t - \boldsymbol{M}_{*} \right)\boldsymbol{V}^{- 1}\left( \boldsymbol{Y}_t - \boldsymbol{M}_{*} \right)^{'}$, $\Tilde{\boldsymbol{B}}=\left( {\alpha_{t}\boldsymbol{\Sigma}}_{L} + \boldsymbol{\Sigma}_{*} \right)$ and $\boldsymbol{\Upsilon}(\alpha_t)=\left( \boldsymbol{\Sigma}_{L} + \boldsymbol{\Sigma}_{*} \right)^{- 1}\boldsymbol{\Sigma}_{*}\Tilde{\boldsymbol{B}}^{- 1}$. \end{enumerate}

Note that the second and third elements within the trace in the derivative at iii) are positive. In contrast, the first one is negative and, for each $\boldsymbol{Y}_t$, $\partial_{\alpha_t}H_t(\alpha_t)\rightarrow -\infty$ as $\alpha_{t}\rightarrow 0^{+}$ and $\partial_{\alpha_t} H_t(\alpha_t)\rightarrow (-n\text{tr}\text(\left( \boldsymbol{\Sigma}_{L} + \boldsymbol{\Sigma}_{*} \right)^{- 1})+\text{tr}(\boldsymbol{\Upsilon}(1)\Tilde{\boldsymbol{A}}))/2$ as $\alpha_{t}\rightarrow 1^{-}$. Thus there is a change of the sign, provided $\text{tr}(\left( \boldsymbol{\Sigma}_{L} + \boldsymbol{\Sigma}_{*} \right)^{- 1}\boldsymbol{\Sigma}_{\ast}(\left( \boldsymbol{\Sigma}_{L} + \boldsymbol{\Sigma}_{*} \right)^{- 1}\tilde{A}-n \boldsymbol{I}))>0$, and the BF has at least one stationary point.

proposition[BF properties, unknown $\boldsymbol{V}$] Under Assumptions (ref) and (ref), the BF $H_{t}(\alpha_t)$ satisfies the following properties. \begin{enumerate}[i)] • There exists a function $\kappa_t\left( \alpha_{t} \right)= (k_{*}k_{A,d}(k_{d}k_{A,*}))^{-np/2}$ $\Gamma_{n}((m_{d}-n-1)/2)$ $\Gamma_{n}((m_{A,*}-n-1)/2)$ $(\Gamma_{n}((m_{*}-n-1)/2)\Gamma_{n}((m_{A,d}-n-1)/2))^{-1}$ which does not depend on $\boldsymbol{Y}_t$ such that $H_t(\alpha_t)\leq \kappa_t(\alpha_t)$, $\kappa_t$ is decreasing in $\alpha_t$, $\kappa_t \rightarrow 1$ as $\alpha_{t} \rightarrow 1^{-}$, $\kappa_t \rightarrow +\infty$ as $\alpha_{t} \rightarrow \underline{\alpha}^{+}$ where $\underline{\alpha}=(2n+p)/m_d$; • $H_{t}\left( \alpha_{t} \right) \rightarrow 1$ as $\alpha_{t} \rightarrow 1^{-}$ and $H_{t}\left( \alpha_{t} \right)\rightarrow\infty$ as $\alpha_{t} \rightarrow \underline{\alpha}^{+}$; • $ \partial_{\alpha_t}H_t(\alpha_t)=L(|\boldsymbol{\Psi}_{A,*}|^{(m_{A,*}-n-1)/2}k_{A,*}^\frac{np}{2}\Gamma_n((m_{A,d}-n-1)/2))^{-2}(B_1 \partial_{\alpha_t}A_1 -A_1 \partial_{\alpha_t}B_1)$, where $L=\left|\boldsymbol{\Psi}\right|_*^{(m_*-n-1)/2}k_*^{np/2}\Gamma_n\left((m_d-n-1)/2)\right)(k_d^{np/2}\left|\boldsymbol{\Psi}_d\right|^{(m_d-n-1)/2}\Gamma_n\left((m_*-n-1)/2\right))^{-1}$, $A_1=a_1a_2a_3$, $B_1=b_1b_2b_3$ with $a_1=k_{A,d}^{np/2}$, $a_2=\Gamma_n\left((m_{A,*}-n-1)/2\right)$ $a_3=\left|\boldsymbol{\Psi}_{A,d}\right|^{(m_{A,d}-n-1)/2}$, $b_1=k_{A,*}^{np/2}$, $b_2=\Gamma_n\left((m_{A,d}-n-1)/2\right)$, $b_3=\left|\boldsymbol{\Psi}_{A,*}\right|^{(m_{A,*}-n-1)/2}$. \end{enumerate}

In the derivative $\partial_{\alpha_t}H_t(\alpha_t)$, the terms $a_i$, $b_i$, $i=1,\ldots,3$ are positive, and the terms $\partial_{\alpha_t}A_1$ and $\partial_{\alpha_t}B_1$ can take both positive and negative values since they are sums of digamma functions. Thus, for some values of $\alpha_t$, the $\partial_{\alpha_t}H_t(\alpha_t)$ can go to zero, and the BF has at least one stationary point.

The following remark discusses the relationship with the univariate case, whereas the empirical section will provide further illustrations for the matrix case.

rem[Univariate Gaussian model] The BF for the univariate outlier detection and its upper bound $\kappa_t(\alpha_t)$ given in Remark (ref) can be easily obtained setting $n = p = 1$ and $\boldsymbol{\Sigma}_{P} = \sigma^{2}/\varphi$ in Eq. (ref), in Prop. (ref) and in i) of Prop. (ref). Since $A_t$ in Remark (ref) is strictly negative, then the upper bound satisfies the properties in Prop. (ref). The derivative of the BF is \begin{equation} \partial_{\alpha_t} H_t(\alpha_t)=\frac{H_{t}(\alpha_t)}{2}\left(\frac{(\varphi+t-1)(Y_t-m_{\ast})^2}{\sigma^2 (\alpha_t(\varphi+t-1)+1)^2}-\kappa_t(\alpha_t)\frac{1}{\alpha_t^2(\varphi+t)}\right) \end{equation} which goes to $-\infty$ as $\alpha_t\rightarrow 0^{+}$ and is strictly positive for $\alpha_t>1/((\varphi+t-1)((Y_t-m_{\ast})^2/\sigma^2-1))$ provided that $(Y_t-m_{\ast})^2/\sigma^2>(\varphi+t)/(\varphi+t-1)$. This implies there exists a stationary point and that $H_{t}<1$ for some values of $\alpha_t$ provided the conditions in Prop. (ref) are satisfied.

Let us now find the subset of the sample space such that for a given threshold $0 < h_{0}(\alpha_t) < \kappa_t(\alpha_t)$, the BF leads us to accept the null hypothesis, that is, $H_{t} \geq h_{0}(\alpha_t) $. Define $\mathbf{y}_{t} = \text{vec}\left( \boldsymbol{Y}_t \right)$ and $\mathbf{m}_{*} = \text{vec}\left( \boldsymbol{M}_{*} \right)$ and let $\kappa_t(\alpha_t)$ be the upper bound given in Prop. (ref). From the expression of the BF given in Prop. (ref) the condition $H_{t} \geq h_{0}(\alpha_t)$ is satisfied for $\mathbf{y}_t$ in the ellipsoid:

equation[equation omitted — 267 chars of source]

centered in $\mathbf{m}_{*}$, with axis in the direction of the eigenvector $\boldsymbol{\nu}_{k}$ of $\boldsymbol{\Sigma}_{H} \otimes \boldsymbol{V}$, where $\boldsymbol{\Sigma}_H=(1 - \alpha_{t})^{-1}\left( \alpha_{t}\boldsymbol{\Sigma}_{L} + \boldsymbol{\Sigma}_{*} \right)\boldsymbol{\Sigma}_{*}^{- 1}\left( \boldsymbol{\Sigma}_{L} + \boldsymbol{\Sigma}_{*} \right)$. The length of the ellipsoid axes along the $\nu_{k}$ eigenvector's direction is $\ell_{k} = \ 2\sqrt{2\log\left( {\frac{\kappa_t(\alpha_t)}{h_{0}(\alpha_t)}} \right)\xi_{k}}$, where $\xi_{k}$ is the corresponding eigenvalue of $\boldsymbol{\Sigma}_{H} \otimes \boldsymbol{V}$. Let $\gamma_{i}$ and $\boldsymbol{\zeta}_{i}$ be the eigenvalues and the corresponding eigenvectors of $\boldsymbol{\Sigma}_{H}$ and $\tau_{j}$ and $\boldsymbol{\delta}_{j}$ the eigenvalues and eigenvectors of $\boldsymbol{V}$. Then $\boldsymbol{\Sigma}_{H} \otimes \boldsymbol{V}$ has eigenvalues $\xi_{k} = \gamma_{i}\tau_{j}$ with corresponding eigenvectors $\boldsymbol{\nu}_{k}= \boldsymbol{\zeta}_i \otimes \boldsymbol{\delta}_{j}$. If we consider $h_{0} = 1$, then $\boldsymbol{Y}_t$ is not considered an outlier for $H_{t} \geq 1$. Figure (ref) provides a numerical illustration for $p=2$ and $n=1$. Increasing the discounting parameter value reduces the evidence in favor of the null for $h_0>1$ (dashed lines, left plot) and increases it for $h_0<1$ (solid lines). For larger values of $h_0$, the evidence against the null becomes stronger (right plot). This effect can also be understood from the following remark on the univariate outlier detection, a special case of Eq. (ref).

figure[figure omitted — 1,010 chars of source]

rem[Univariate Gaussian model] In the univariate case, the null hypothesis is accepted, i.e. $H_t>h_{0}$, for values of $Y_t$ in the interval $(Y_{1t},Y_{2t})$ with \begin{equation} Y_{jt}=m_{\ast}+(-1)^{j}\sqrt{\log\left(\frac{h_{0}(\alpha_t)}{\kappa_t(\alpha_t)}\right)\frac{2\sigma^2(\varphi+t)(\alpha_t(\varphi+t-1)+1)}{(\alpha_t-1)(\varphi+t-1)}}. \end{equation}

The BF and the outcome of the testing procedure depend on the choice of the discounting parameter $\alpha_t$. In this paper, we propose the integrated BF and its normalized version as a solution, which requires that the integral with respect to $\pi_t(\alpha_t)$ is bounded. In the following, we provide existence conditions for the IBF and NIBF under the beta perturbation assumption, a standard prior distribution used in Bayesian inference for parameters on bounded intervals. As stated in the following proposition, the integral of the upper bound for the univariate case can be derived analytically and is well-defined under the assumption of a standard uniform distribution for the discounting parameter.

proposition[Integrability] Assume a general beta distribution $\mathcal{B}e(a,b)$ truncated on the interval ($\underline{\alpha},\bar{\alpha})$. The IBF, $H_t=\int_{0}^{1}H_t(\alpha)\pi_{t}(\alpha)d\alpha$, satisfies the following properties. \begin{enumerate} • Under Assumptions (ref) and (ref), for $a>np/2$, $b>0$, $\underline{\alpha}=0$ and $\bar{\alpha}=1$ \begin{equation} H_t<\sum_{\ell=0}^{\infty}\sum_{k\in\mathcal{K}_{\ell}}d_{\ell,k}\bar{c}_{k}\frac{\Gamma(n/2+1)B(a-np/2+w_{k},b)}{\Gamma(\ell+1)\Gamma(n/2-\ell)B(a,b)}<\infty \end{equation} with $\underline{k}=(k_0,\ldots,k_p)$, $\mathcal{K}_{\ell}=\{\underline{k}|k_0+\ldots,k_p=\ell\}$, $\bar{c}_{0}=\tilde{c}_{0}-1$, $\bar{c}_{j}=\tilde{c}_j$, $j=1,\ldots,p$, and $$ d_{\ell,\underline{k}}=\binom{\ell}{k_0,k_1,k_2,\ldots,k_p},\quad \bar{c}_{\underline{k}}=\prod_{j=0}^{p}\bar{c}_{j}^{k_j}, $$ where $d_{\ell,\underline{k}}$ denotes the multinomial coefficient, $\tilde{c}_j$ $j=0,\ldots,p$ are the coefficients of the characteristic polynomial of $\boldsymbol{\Sigma}_L^{-1}\boldsymbol{\Sigma}_{\ast}$. • Under Assumptions (ref) and (ref), for $a>1$ and $b>0$, $\underline{\alpha}=(2n+p)/m_d$ and $\bar{\alpha}=1$. \begin{align} &H_t<(\alpha k_{*}/(\alpha(k_{*}+1)))^{np/2}\Gamma\left(\frac{p}{2}\right)^{-n}\left(\frac{4}{p}\right)^n\frac{B(1-\underline{\alpha};ra-r+1,qb-q+1)}{B(1-\underline{\alpha};a,b)}\xi(r,q) \end{align} where $B(c,a,b)=\int_{0}^{c}x^{a-1}(1-x)^{b-1}dx$ denotes the incomplete beta function, $\xi(r,q)=(((m_{d}(1-\underline{\alpha})+p)^{1-qn}-p^{1-qn})/(m_{d}(1-qn)))^{1/q}$ and $q=r/(r-1)$ with $r > 1$. \end{enumerate}

corollary[IBF in the univariate model] The integral of the BF $H_t$ given in Remark (ref) is bounded and satisfies: \begin{equation} \int_{0}^{1}{H_{t}\left( \alpha_{t} \right)d\alpha_{t} - 1 \leq \frac{1}{2\sqrt{(\varphi + t-1)(\varphi + t)}}\log\left( \frac{\sqrt{(\varphi + t)} + \sqrt{(\varphi + t-1)}}{\sqrt{(\varphi + t)} - \sqrt{(\varphi + t-1)}} \right)}. \end{equation}

In the empirical illustration, we show that the beta hyper--parameters' choice can affect the hypothesis testing outcome, which calls for calibrated BFs. We assume $\boldsymbol{Y}_t$ in the predictive BF is not observed at time $t$. Thus, the predictive BF is random, and we prove that its distribution is a mixture of gamma distributions. Controlling for the test's size and power while using the BF distribution allows us to derive a calibrated BF and a suitable critical region for the test.

proposition[Distribution of the BF, known $\boldsymbol{V}$] Assume $\boldsymbol{Y}_t|\mathcal{D}_{t-1}\sim\mathcal{N}_{p,n}(\tilde{\boldsymbol{M}},\tilde{\boldsymbol{\Sigma}},\boldsymbol{V})$, let $\lambda_1$, $\ldots,\lambda_p$ be the eigenvalues of the covariance matrix $(\boldsymbol{\Sigma}_{H}^{-1/2})'\tilde{\boldsymbol{\Sigma}}\boldsymbol{\Sigma}_{H}^{-1/2}$ with $\boldsymbol{\Sigma}_H= (1 - \alpha_{t})^{-1}\left( \alpha_{t}\tilde{\boldsymbol{\Sigma}} + \boldsymbol{\Sigma}_{*} \right)\boldsymbol{\Sigma}_{*}^{- 1}\left( \tilde{\boldsymbol{\Sigma}} + \boldsymbol{\Sigma}_{*} \right)$ and $Q$ the associated eigenvector matrix. Under Assumptions (ref) and (ref), the distribution of the BF $H_t$ defined in Eq. (ref) has pdf and cdf: \begin{eqnarray} f_{H_{t}}(h|\alpha_t) &=&\frac{2}{h}\sum_{k=0}^{\infty}c_k g(- 2\log\left( h/\kappa_{t} \right);\frac{np}{2}+k,2\lambda),\\ F_{H_{t}}(h|\alpha_t) &=&\sum_{k=0}^{\infty}c_k (1-G(- 2\log\left( h/\kappa_{t} \right);\frac{np}{2}+k,2\lambda)), \end{eqnarray} respectively, with support $0<h\leq\kappa_{t}$, where $\kappa_t(\alpha_t)$ is defined in Prop. (ref), the parameter $0<\lambda<\infty$ is arbitrarily chosen, $g(x;a,b)$ and $G(x;a,b)$ are the pdf and cdf of a gamma distribution with shape and scale parameters $a>0$ and $b>0$, respectively, and the coefficients $c_k$, $d_k$ and $f_k$ satisfy \begin{eqnarray*} c_{k}&=&\exp(-\sum_{j=1}^p U_{jj})\prod_{j=1}^{p} (\lambda_j/\lambda)^{-n/2}f_k,\quad d_k=\frac{n}{2k}\sum_{j=1}^{p}(1-\lambda/\lambda_j)^{k} + \lambda \sum_{j=1}^{p}\frac{U_{jj}}{\lambda_j}(1-\lambda/\lambda_j)^{k-1},\\ f_{k+1}&=&\frac{1}{k+1}\sum_{j=1}^{k+1}j d_j f_{k+1-j}, k=0,1,2,\ldots,\quad with\quad f_0=1, \end{eqnarray*} where $U_{jj}$ are the $(j,j)$-th entry of $\boldsymbol{U}=\boldsymbol{Q}'(\boldsymbol{\Sigma}_{H}^{1/2}\tilde{\boldsymbol{\Sigma}}^{-1}(\tilde{\boldsymbol{M}}-M_{\ast})\boldsymbol{V}^{-1})(\tilde{\boldsymbol{M}}-\boldsymbol{M}_{\ast})'\boldsymbol{\Sigma}^{-1/2}_{H}\boldsymbol{Q}$.

If in the previous proposition, we set $\tilde{\boldsymbol{M}}= \boldsymbol{M}_{*}$ and $\tilde{\boldsymbol{\Sigma}}= \boldsymbol{\Sigma}_{d}$, that is the location and scale of the marginal likelihood under the null (from i) in Prop. (ref)) we denote the distribution with $F_{0,t}(h)$. In contrast, if $\tilde{\boldsymbol{M}}= \boldsymbol{M}_{*}$ and $\tilde{\boldsymbol{\Sigma}}= \boldsymbol{\Sigma}_{A,d}$, that is the location and scale of the marginal likelihood under the alternative (from i) in Prop. (ref)), we denote the distribution with $F_{1,t}(h)$.

From the properties of the gamma distribution, $f_{H_{t}}(h|\alpha_t)$ is continuous at $0$ and $\kappa_t$. As stated in the following, the result of Prop. (ref) provides the BF distribution for the univariate Gaussian model given in previous studies weiss1997bayesian,DESANTIS2004121,pawel2025closed.

corollary[Univariate Gaussian model] In the univariate case, i.e. $n=p=1$, following the notation in Remark (ref), we set $\boldsymbol{\Sigma}=\sigma^2$, $\boldsymbol{V}=1$ and $\boldsymbol{B}=\theta$. Assuming $\tilde{\boldsymbol{M}}=\boldsymbol{B}$ and $\tilde{\boldsymbol{\Sigma}}=\boldsymbol{\Sigma}_L$ in Prop. (ref) the distribution of the BF is \begin{equation} F_{H_{t}}\left( h|\alpha_{t} \right) =1-\Phi(\sqrt{-2\log(h/\kappa_t)\sigma^2_H/\sigma^2_L}-\sqrt{\gamma})+\Phi(-\sqrt{-2\log(h/\kappa_t)\sigma^2_H/\sigma^2_L}-\sqrt{\gamma}), \end{equation} where $\gamma=(\theta-m_{\ast})^2/\sigma^2$, the variances $\sigma^2_H=\boldsymbol{\Sigma}_H$ and $\sigma^2_L=\boldsymbol{\Sigma}_L$ are given in Prop. (ref), the bound $\kappa_t$ is given in Remark (ref) and $\Phi(\cdot)$ is the cdf of a standard normal distribution.

Our perturbation framework for matrix--variate observations is a form of global regularization that affects all entries within the matrix. Later in this paper, we will explore the sensitivity of the testing procedure to variations in the patterns, proportions, and magnitude of outliers. As an extension, local contamination frameworks can be devised to detect patterns within the outliers and potential dependencies among outlying observations. In the context of matrices, a multiplicative perturbation can be employed to define $p_A(\boldsymbol{Y}_t|\mathcal{D}_t)$, assuming, for instance, the perturbed normal distribution $\mathcal{N}_{p,n}(M_{\ast}, \boldsymbol{A}_1\boldsymbol{\Sigma}_d \boldsymbol{A}_1, \boldsymbol{A}_2 \boldsymbol{V} \boldsymbol{A}_2)$, where $\boldsymbol{A}_1$ and $\boldsymbol{A}_2$ are two diagonal matrices scaling along the rows and columns with varying levels. An additive perturbation would lead to a distribution like $\mathcal{N}_{p,n}(\boldsymbol{M}_{\ast},\boldsymbol{\Sigma}_d+\boldsymbol{A}_1,\boldsymbol{V}+\boldsymbol{A}_2)$, with $\boldsymbol{A}_1$ and $\boldsymbol{A}_2$ shifting the posterior's covariance matrix diagonal. On the one hand, local contamination allows for patterns among outliers; however, it requires the identification of these patterns, which in turn calls for inference on the contamination parameters. Therefore, alternative approaches, such as multiple-outlier models presented in page2011bayesian,tomarchio2022mixtures, might be preferable. This typically involves introducing appropriate prior specifications and has the drawback of requiring computationally intensive procedures for posterior approximation. We will postpone this extension to future research.

figure[figure omitted — 1,978 chars of source]

Simulation study

In the following simulation exercise, we generate data with outliers and study some finite-sample properties of the decision procedures based on standard BF, MBF, IBF and NIBF, and the calibrated predictive BF approach we introduced. All analyses were implemented in Matlab 2023a and carried out on a 12-core computing system with 128GB RAM. The parameter values have been randomly generated as follows: $\boldsymbol{M}\sim \mathcal{N}_{p,n}(\boldsymbol{O}_{p\times n}, \boldsymbol{I}_{p},\boldsymbol{I}_{n})$, $\boldsymbol{\Sigma}=\boldsymbol{SS}'$, $\boldsymbol{S}\sim\mathcal{N}_{p,p}(\boldsymbol{O}_{p\times p}, \boldsymbol{I}_{p},\boldsymbol{I}_{p})$, $\boldsymbol{\Psi}=\boldsymbol{GG}'$, $G\sim\mathcal{N}_{n,n}(\boldsymbol{O}_{n\times n}, \boldsymbol{I}_{n}, \boldsymbol{I}_{n})$. We consider two scenarios: the case of moderate-size observation matrix, denoted as $Case_1$, with $p = 30$ and $n = 10$, and the case of a large-size observation matrix, $Case_2$, with $p = 50$ and $n = 50$. We generate the synthetic data with an outlier as follows. First a noise sequence is generated, $\boldsymbol{E}_{t}\sim\mathcal{N}_{p,n}(\boldsymbol{O}_{p\times p},\boldsymbol{\Sigma},\boldsymbol{\Psi})$ i.i.d. $t=1,\ldots,100$. Secondly, the observable sequence is defined as $\boldsymbol{X}_{t}=\boldsymbol{M}+\boldsymbol{E}_t$ for $t\neq 80$ and $\boldsymbol{X}_{t}=\boldsymbol{M}+u \boldsymbol{R}_t+\boldsymbol{E}_t$ for $t=80$, where $\boldsymbol{R}_t$ is a binary matrix encoding the position of the outliers and $u$ is 0.5, 1, 1.5, 3, 5 and 15 that are $1/30$, $1/15$, $1/10$, $1/5$, $1/3$ and $1$ average standard deviation of the matrix normal.

Our simulation study considered various settings obtained by varying the configurations of the matrix $\boldsymbol{R}_t$ and the magnitude of the outliers. For each setting $J$ independent datasets have been generated $\boldsymbol{X}_{t}^{(j)}$, $t=1,\ldots,T$ for $j=1,\ldots,J$, with $T=100$ and the BFs $H_t^{(j)}$ computed. The probabilities $p_{I}=P(H_t>\bar{h})$, $p_{II}=P(H_t<\underline{h})$ and $p_{III}=P(\bar{h}<H_t<\underline{h})$ have been estimated as follows

equation*[equation* omitted — 259 chars of source]

under the null hypothesis of the absence of outliers at $t\neq 80$ and under the alternative hypothesis of a certain number of outlying observations in the observation matrix at $t=80$. A summary of the results for $Case_1$ and $Case_2$ is provided in Tab. (ref) and Tab. (ref) of Appendix (ref) respectively. The results are obtained with $J=100$ experiments for the two outlier settings. The two settings differ for the positions of the outliers within the matrix $R_t$: i) row and column patterns in the positions, panel (a); and ii) completely random positions, panel (b). Appendix (ref) provides reference to an R package we specifically developed for outlier detection.

Figure (ref) illustrates the result of the BF procedure on one of the simulated datasets for $u=0.5$ (top left panel a in Tab. (ref)). In the BF procedure the calibrated value of $\alpha_t$ for $\tau=0.01$ and $\beta=0.8$ is $\alpha^{*}=0.750$ for all $t$ and the inconclusive interval $C_{H}$ has lower and upper bounds $\underline{h}(\alpha^{*})=0.839$ $\bar{h}(\alpha^*)=1.161$, respectively. Starting from the left in the first row, we notice that the observations with a change in the convexity of $H_t(\alpha_t)$ are recognized as outliers based on the value of $\alpha_t$. Notably, for the outlier observation introduced in the simulation, the BF is consistently convex and remains far below 1. The second plot illustrates the behavior of the thresholds $\underline{h}(\alpha_t)$ and $\bar{h}(\alpha_t)$ defined at the end of Section (ref) for given power and size values. The length of the randomized decision interval $C_h$ (vertical axis) reduces as $\alpha_t$ goes to 1. The third plot presents the BF and the inconclusive region for $\alpha^{*}=0.750$. The BF of the outlying observation is located outside the inconclusive interval, whereas about 16 observations exhibit a BF below 1 within the inconclusive interval. Moving to the second row, the Minimum BF suggests positive evidence against the null hypothesis of the absence of outliers ($MBF<10^{-1/2}$ following Jeffrey's scale of evidence) for three observations. The IBF and NIBF (solid line in the second and third plot at the bottom) provide strong evidence against the null ($IBF$, $NIBF<10^{-1/2}-1\approx -0.687$) only for the 80th observation and barely worth mentioning evidence ($IBF$, $NIBF>10^{-1/2}-1\approx -0.687$) for the other observations. Nevertheless, comparing the solid and dashed lines shows that the outcome of the IBF and NIBF procedures strongly depends on the choice of the hyperparameters of the beta distribution. Thus, the third robust method should be applied, exploiting the BF's sampling variability to find the calibrated $\alpha_t$ and the reference thresholds.

In all the experiments, when data are generated under the null, the type I error probability $P(H_t<\underline{h}|\mathcal{H}_0)$ is about 2%, whereas when data are generated under the alternative hypothesis (first column), the power $P(H_t<\underline{h}|\mathcal{H}_1)$ of the test gets close to one, increasing the number of outliers (e.g., see columns of the panels (b) in Tab. (ref)). As we can expect, the convergence is faster for larger number of observations (see panels (b) in Tab. (ref)). Also, the effective size and power may depend on the position of outliers in the rows or columns of the matrix $R_t$. The presence of patterns in the position of the outliers reduces the power compared to the case of a completely random position within the matrix $R_t$ (compare columns $20\times 10$ and $200$ in panels (a) and (b), respectively). The power decreases below 80% for small outlier amplitude (e.g. $1/10$ standard deviation) and a small number of outlying observations (e.g. 10 out of 300 elements) within the matrix. When all entries are outliers, the power of the test converges to 1, increasing the outliers' magnitude (different rows in the last column in Tab. (ref)).

Empirical Illustration

We illustrate our sequential matrix outlier detection on three relevant benchmark datasets. See Appendix (ref) for a detailed description of the datasets and a discussion of the outliers' dating.

Inflation and Unemployment Dataset

The Inflation and Unemployment Dataset (Can09) spans a period from January 2002 to October 2022 ($T =250$) and consists of a sequence of $p\times n$ matrices, covering $p = 11$ EU countries and $n=3$ macroeconomic variables (the Industrial Production Index, the Price Index and the Unemployment Rate,). The total number of observations is $8,250$, representing an example of big data in this field.

In the testing procedure with BF, we consider rolling windows of $w$ observations each and assumed $\boldsymbol{\Sigma}_P=\boldsymbol{\Sigma}_L /\varphi$, with $\varphi=w$. For each window, the prior mean $\boldsymbol{M}$ is set equal to the posterior mean $M_{\ast}$ of the previous window. The variance-covariance $\boldsymbol{\Sigma}_L$ has been estimated by Least Squares.

The top plots in Figure (ref) display the BF $H_t(\alpha_t)$ (solid grey) and the upper bound of the BF $\kappa_t(\alpha_t)$ (solid red) as functions of $\alpha_t$ for various dates (each represented by a different line). The findings in the top-left plot suggest that it is crucial to compare the BF to the upper bound. The BF tends to hover around one for higher values of $\alpha_t$, with some exceptions, and it significantly exceeds one only for small values of $\alpha_t$ (illustrated by the dark grey lines). It frequently intersects the threshold at lower values of $\alpha_t$ (light grey) as well as for higher values. Additionally, there are instances where it crosses the threshold intermittently.

The left-bottom plot offers a different illustration of the effect of $\alpha_t$ on the outcome of the sequential outlier detection for the entire sample from September 2018 to October 2022.

figure[figure omitted — 1,885 chars of source]

In all settings, an outlier was detected in March 2020 (i.e. at the pandemic outbreak). For $\alpha_t=0.054$ (solid line), a sequence of outliers is detected before February 2020 with BF far from 1, whereas, for $\alpha_t$ equal to 0.402, 0.801 and 0.851 (dashed, dotted and dashed-dotted, respectively), the BF is close to one before March. Qualitatively speaking, after the outbreak, for $\alpha_t=0.402$, the BF is close to one after December 2021, whereas larger $\alpha_t$ values return BF close to just one after March 2020. For illustrative purposes, we present in detail the results for some relevant dates (Figure (ref) in the Appendix). For the observations in February 2020, during the pandemic outbreak, the null hypothesis of the absence of an outlier is not rejected for any choice of $\alpha_t$. Nevertheless, the BF is close to one for large values of $\alpha_t$. In contrast, in March 2020, there is strong evidence of an outlier. On the other dates in the figure, the alternative hypothesis is accepted for some values of $\alpha_t$. The dashed vertical lines indicate the stationary point. The stationary point is not near zero on some dates, such as February 2022. On other dates, the stationary point is near zero, and the BF is far below one for a large part of the $\alpha_t$ values (e.g., March 2020, March 2022 and October 2022). The top-left plot in Figure (ref) shows the results of the alternative procedures.

The MBF and IBFs provide clearer identification of outlying observations, enhancing the standard BF and supporting the outcome of the testing procedure endowed with the inconclusive interval given in the first line of Figure (ref). Dashed lines indicate the lower and upper bounds of the inconclusive region, the dots show the value of the BF. For $\alpha_t=0.75$ and a test size of 1%, the power is approximately 96%.

The outcome of the sequential test is in line with those of classical frequentist tests for outliers such as the Grubb's (G) test Grubb50 and the Generalized Extreme Studentized Deviate (GESD) test rosner1983percentage. The second row of Figure (ref) reports the number of outliers detected at the 1% level by applying the two tests element-wise to each entry of the observation matrices.

figure[figure omitted — 2,591 chars of source]

figure[figure omitted — 1,685 chars of source]

Trade Network Dataset

We consider the Trade Network Dataset (Rose2004), provided by the IMF, as it is a key reference in international trade studies. It integrates country reports with data from COMTRADE and EUROSTAT, covering 159 countries from 1995 to 2017 with annual frequency. The dataset sample we consider includes a sequence of 22 import networks of trade across 27 countries, i.e. $n=p=27$.

Following the results in the top-middle plot of Figure (ref), we find evidence of outlier observations on all dates except for 2016. Nevertheless, the BF is close to one in some cases, such as 2015 and 2017. From the sensitivity analysis in Figure (ref), one can see that the test outcome depends crucially on the choice of $\alpha_t$. The top-middle plot in Figure (ref) indicates that minimum and integrated BFs provide evidence of the absence of an outlier in 2016 and the presence of outlying observations in 2015 and 2017, supporting the conclusion of the calibrated BF procedure with inconclusive intervals presented in the Panel (a) middle plot of Figure (ref). In the year 2016, the classical G and GESD tests for outliers, applied entry-wise to the observation matrix, detected one outlying series out of 729 series (mid plot in Panel (b).

Volatility Network Dataset

We investigate the presence of outliers in a sequence of volatility networks among European firms with the largest market capitalization (billio2021matrix). The dataset consists of 145 temporal networks (from the 4th of January 2016 to the 30th of September 2020) between 50 firms, i.e. $T=145$ and $n=p=50$, for 362,500 observations. In the sequential outlier detection, we used a rolling window of 90 observations.

The plots in the top-right and bottom-right sections of Figure (ref) indicate that the rejection of the hypothesis is independent of the chosen discounting coefficient. However, BF exhibits greater sensitivity to discounting towards the end of 2019. The sensitivity analysis of the BF to $\alpha_t$ is presented for selected dates in Figure (ref) in the Appendix.

The top-right plot in Figure (ref) shows the results of the minimum BF, which support the main findings of the calibrated BF procedure (Figure (ref), Panel(a), right plot). When the BF is far above one, the classical G and GESD tests for outliers, applied entry-wise to the observation matrix, detect a reduction in the number of outliers. In summary, the rapid changes in volatility and the persistence of volatility regimes call for nonlinear models, such as switching or threshold models, which account for structural breaks and recurrent regimes.

Conclusion

The assessment of the model performance is relevant in many applications and becomes crucial in forecasting. This paper proposes sequential outlier detection for the matrix normal model. The hypothesis testing procedure extends the predictive Bayes Factor (BF) with the power discounting to matrix models. The proposed approach relies on normality, now a default assumption in many applications, which serves to build a preliminary test for outliers in sequences of matrix--valued data. Some solutions are proposed to mitigate the test outcome's dependence on the discounting coefficient value, such as the minimum and the integrated BFs. The finite--sample distribution of the predictive BF is derived, and a testing procedure is proposed based on calibrated discounting and BF. Simulation experiments are conducted to study the properties of our Bayesian outlier detection. Numerical illustrations with relevant benchmark datasets are given. They include a comparison with classical tests for outliers and a validation based on major global event dates.

\phantomsection{ Funding}

This work was funded by the MUR -- PRIN project under g.a. n. 2022CLTYP4 and the Next Generation EU -- `GRINS -- Growing Resilient, INclusive and Sustainable' project (PE0000018), National Recovery and Resilience Plan -- PE9. The views and opinions expressed are only those of the authors and do not necessarily reflect those of the EU.

\phantomsection

{\bf Supplementary Materials}

description• A PDF document containing all mathematical derivations, proofs of theoretical results, and additional simulation evidence supporting the methodology presented in the main text. (PDF file) • The R package BAYSFWATCH, developed for this study, contains functions to implement the proposed Bayesian outlier detection procedures. The package also includes all datasets used in the empirical illustrations. (GNU zipped tar file, available at \url{https://github.com/BayesianEcon/BAYSFWATCH})

center[center omitted — 110 chars of source]