# Limitations of Out-of-Distribution Detection in 3D Medical Image
Segmentation

Anton Vasiliuk, Daria Frolova, Mikhail Belyaev, and Boris Shirokikh The
results have been obtained under the support of the Russian Science Foundation
grant number 20-71-10134. (Corresponding author: B. Shirokikh.)A. Vasiliuk is
with the Moscow Institute of Physics and Technology, Dolgoprudny, Moscow
Region, 141701, Russia and Artificial Intelligence Research Institute (AIRI),
Moscow, 105064, Russia (e-mail: vasilyuk@phystech.edu).D. Frolova, M. Belyaev,
and B. Shirokikh are with the Skolkovo Institute of Science and Technology,
Moscow, 143026, Russia and Artificial Intelligence Research Institute (AIRI),
Moscow, 105064, Russia (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>boris.shirokikh@skoltech.ru).This work has been submitted to the IEEE for
possible publication. Copyright may be transferred without notice, after which
this version may no longer be accessible.

###### Abstract

Deep Learning models perform unreliably when the data comes from a
distribution different from the training one. In critical applications such as
medical imaging, out-of-distribution (OOD) detection methods help to identify
such data samples, preventing erroneous predictions. In this paper, we further
investigate the OOD detection effectiveness when applied to 3D medical image
segmentation. We design several OOD challenges representing clinically
occurring cases and show that none of these methods achieve acceptable
performance. Methods not dedicated to segmentation severely fail to perform in
the designed setups; their best mean false positive rate at 95% true positive
rate (FPR) is 0.59. Segmentation-dedicated ones still achieve suboptimal
performance, with the best mean FPR of 0.31 (lower is better). To indicate
this suboptimality, we develop a simple method called Intensity Histogram
Features (IHF), which performs comparable or better in the same challenges,
with a mean FPR of 0.25. Our findings highlight the limitations of the
existing OOD detection methods on 3D medical images and present a promising
avenue for improving them. To facilitate research in this area, we release the
designed challenges as a publicly available benchmark and formulate practical
criteria to test the OOD detection generalization beyond the suggested
benchmark. We also propose IHF as a solid baseline to contest the emerging
methods.

###### Index Terms:

Computed Tomography, Magnetic Resonance Imaging, out-of-distribution
detection, segmentation.

##  I Introduction

In recent years, Deep Learning (DL) methods have achieved human-level
performance in automated medical image processing. But the development of
these methods on a large scale is slowed by several factors. One such factor
is the unreliable performance of DL models when the data comes from a
distribution different from the training one [1]. These differences are common
in medical imaging: population, demographic, acquisition parameter changes, or
new imaging modalities.

Out-of-distribution (OOD) detection helps to identify the data samples with
such differences, hence increasing the reliability and safety of a DL model.
For instance, detected cases could be marked as rejected, preserving the model
performance, or reported to the experts, preventing the model from failing
silently. The ability to report or reject unreliable cases is now considered a
necessary capability to enable safe clinical deployment [2].

OOD detection on natural images is a well-researched area [3] where several
established benchmarks [4, 5] facilitate its development. Moreover, these
methods directly scale on 2D medical images, resulting in multiple algorithms
[6, 7, 8], and also a benchmark [9]. At the same time, OOD detection on 3D
medical images remains poorly explored, although 3D medical image segmentation
is one of the most addressed tasks in medical imaging [10] with outstanding
practical usefulness, e.g., quantifying anatomical structures, pathologies, or
important biomarkers.

The primary cause of this poor exploration is the lack of datasets and
benchmarks with a correct problem design. For example, one party uses private
data [11], while the other simulates synthetic anomalies that are unlikely to
occur in clinical settings [12]. A study can be limited to a single
distribution shift, e.g., changes in the scanning location [11], lacking the
diversity of setups. Also, a study can be restricted to uncertainty estimation
[13] or anomaly detection [12] methods, leaving the full spectrum of
approaches uncovered. Such issues limit a fair comparison of the proposed
approaches.

In this paper, we investigate the effectiveness of OOD detection when applied
to 3D medical image segmentation, closing the outlined gaps in prior work. To
enable a correct comparison, we thus design a diverse set of challenges using
publicly available data with a downstream segmentation task and simulating
clinically occurring anomaly sources. Besides the problem design, such a study
requires appropriately selected state-of-the-art methods. We note that several
areas, e.g., anomaly detection and uncertainty estimation, share motivation
and methodology with OOD detection. Therefore, we overview all related areas
and, contrary to the previous works, present a complete methodological
coverage.

Extensive evaluation of six selected methods results in our main conclusion:
state-of-the-art OOD detection falls short of achieving optimal performance on
3D medical images. We show that the methods not designed for segmentation
completely fail in most setups, scoring from $0.84$ to $0.59$ False-Positive
Rate (FPR) on average, which is not far below $0.95$ FPR of the random
guessing. (Lower FPR is better.) Two methods specifically designed for 3D
segmentation achieve $0.38$ and $0.31$ mean FPR, further reducing the error
about two times. At the same time, we show that these errors can be reduced
even further with a simple approach.

We show this space for improvement by developing a histogram-based method
called Intensity Histogram Features (IHF). IHF achieves comparable and often
superior results to its competitors with $0.25$ mean FPR. It also scores $0$
FPR in multiple challenges, indicating that the distribution shifts in 3D
medical imaging can often be detected using image intensity histograms, while
the DL-based methods overlook this domain feature. Therefore, we consider
current DL-based OOD detection far from unveiling its full potential and
assume it can be further improved.

Given IHF’s negligible computational costs compared to DL, we suggest it as a
baseline to contest the emerging OOD detection methods. Furthermore, we
propose using the designed challenges as a benchmark for developing new
methods. Correct problem setting, in-depth analysis with simple methods, such
as IHF, and ablation studies on synthetic data confirm that our benchmark
allows to estimate the quality of solving general OOD detection instead of
classifying a priori known anomaly types. Thus, summarizing our contributions,
we outline the following:

  1. 1.

We demonstrate severe limitations of the existing OOD detection methods on 3D
medical images.

  2. 2.

We design and release the corresponding benchmark that can be used as a
starting point for the related research.

  3. 3.

We propose a method, IHF, suggesting it as a solid baseline for OOD detection
on 3D medical images.

Below, we describe the data used in our study and the problem setup (Sec. II).
Then, we review and select state-of-the-art and core methods from the related
fields and also detail IHF (Sec. III). Finally, we present the results (Sec.
V) and discuss the limitations and implications of our study (Sec. VI).

##  II Data

Figure 1: Examples of CT images (representative axial slices) from different
simulated OOD sources in our benchmark. Figure 2: Examples of MRI images
(representative axial slices) from different simulated OOD sources in our
benchmark.

Contrary to the fields of 2D natural and medical images, no established OOD
detection benchmark with a correct problem setting exists for 3D medical
images. For example, [11] used a variety of brain and abdominal CT and MRI
datasets but included private ones. The authors also studied only a single
distribution shift, changes in the scanning location, which does not allow to
estimate the general performance. [12] created an OOD detection challenge,
simulating synthetic anomalies in brain MR and abdominal CT images. However,
their setup lacks a downstream task (e.g., segmentation), so their study is
limited to unsupervised anomaly detection methods. Synthesizing local
corruptions, as in [12], can also lead to evaluation biases, which we show
with our analysis. On the other hand, [13] included datasets with the
segmentation task but limited the considered methods to supervised uncertainty
estimation.

Given the disagreement of setups, their partial problem coverage, or privacy,
we design the OOD detection challenges from scratch following three core
principles:

  * •

We include two large publicly available CT and MRI in-distribution (ID)
datasets to cover the most frequent volumetric modalities.

  * •

We ensure both datasets have a downstream segmentation task, allowing us to
use the full spectrum of methods.

  * •

We select diverse OOD datasets that simulate the clinically occurring sources
of anomalies: changes in acquisition protocol, patient population, or
anatomical region. All these datasets are also publicly available.

We also introduce several medical imaging artifacts as anomaly sources,
synthetically generated as random transformations. Generating synthetic
anomalies is a popular approach, applied to 3D images [12, 13] as well as 2D
ones [4, 5]; this approach also allows us to conduct controlled ablation
studies at different distortion levels.

We made the resulting benchmark publicly
available111https://github.com/francisso/OOD-benchmark and we detail all
related CT and MRI datasets in Sec. II-A and Sec. II-B. The problem setting is
described in Sec. II-C.

###  II-A 3D CT datasets

We construct a total of 6 challenges on CT data, including two synthetic ones.
We give a visual example of data samples in Fig. 1 and detail the ID dataset
and every setup below.

#### ID dataset

As an ID CT dataset, we use LIDC-IDRI [14]. It contains $1018$ chest CT images
with the lung nodules segmentation task. We remove cases without nodules since
they do not contribute to training a segmentation model. Then, we randomly
split the rest $883$ images $4:1$ into the train and test, stratified by the
number of nodules.

####  II-A1 OOD source: scanner

To simulate a covariate shift, we select Cancer500 [15] that has the same
downstream task as the ID dataset but is obtained with different scanners and
acquisition protocols. It contains $979$ chest CT images. We exclude all
images with low resolution (less than $64$ axial slices) and no annotated
nodules, resulting in $841$ images left.

####  II-A2 OOD source: population

To simulate a patient population shift, we use two datasets with the similar
semantic content but different downstream task. These datasets are
Medseg9222https://radiopaedia.org/articles/covid-19-3 and MIDRC [16],
containing 9 and 154 chest CT images, respectively, with COVID-19 cases.
Excluding all non-COVID cases, the merged dataset has 120 images.

####  II-A3 OOD source: location (liver)

To simulate a semantic shift, we select a dataset of the same modality that
focuses on different body region. Here, we use LiTS [17], a dataset with $201$
abdominal CT images.

####  II-A4 OOD source: location (head)

Similarly, we include CT-ICH [18], a dataset with $75$ head CT images.

####  II-A5 OOD source: synthetic (image noise)

We simulate local image corruptions by applying damaging transformations to
the testing cases of ID dataset. Applied to a randomly selected image crop,
these transformations include blurring, changing contrast, or inserting noise.

####  II-A6 OOD source: synthetic (elastic)

We simulate tissue anomalies by applying elastic transform of random severity.

###  II-B 3D MRI datasets

We construct a total of 7 challenges on MRI data, including four synthetic
ones. We give a visual example of data samples in Fig. 2 and detail every
setup below.

#### ID dataset

As an ID MRI dataset, we use VS-Seg [19]. It contains $242$ brain T1c MRIs
with the vestibular schwannoma segmentation task. We remove cases with empty
target mask and split the rest $239$ images $2:1$ into the train and test.

####  II-B1 OOD source: scanner

To simulate a covariate shift, we select the data with the same semantic
content and downstream task but is obtained with different scanners and
acquisition protocols. Here, we choose CrossMoDA ETZ as a subset of the
CrossMoDA 2022 Challenge dataset [20] with $105$ brain T1c MR images and use
it without changes.

####  II-B2 OOD source: population (glioblastoma)

To simulate a patient population shift, we use EGD [21], a dataset with 774
brain MRIs of four modalities (FLAIR, T1, T1c, T2) with a glioma segmentation
task. We reduce a possible covariate shift by using only T1c modality from the
Siemens Avanto 1.5T scanner, as in VS-Seg, resulting in 262 seleted images.

####  II-B3 OOD source: population (healthy)

Additionally, we simulate a patient population shift with healthy cases
instead of the changing pathology. To do so, we use the CC359 [22] dataset
with $359$ brain MR images of T1 modality. We note however that CC359 images
differ in vendor, scanning protocol, and do not contain contrast enhancement,
so this setup has a secondary OOD source, a covariate shift.

####  II-B4 OOD source: synthetic (K-space noise)

We synthesize MR imaging artifact, known as Herringbone artifact, at different
magnitudes. This results in the visible spikes across the whole image due to
anomaly points in K-space.

####  II-B5 OOD source: synthetic (anisotropy)

We synthesize wrong resolution by downsampling the image and upsampling it
back along one randomly chosen axis.

####  II-B6 OOD source: synthetic (motion)

We synthesize two types of MR imaging artifacts that can happen due to the
patient motion. One is ghosting, which appears as shifted copies of the
original image. Another exploits RandomMotion simulation from the torchIO
library [23].

####  II-B7 OOD source: synthetic (image noise)

The same pipeline as for CT images; see Sec. II-A5.

###  II-C Problem setting

We define the OOD detection problem as the classification between samples from
a source distribution (ID) and abnormal samples from a novel different
distribution (OOD). The core assumption is the abnormal sample distribution is
unknown and cannot be computed in advance. Thus, we approximate the anomaly
distribution by constructing a diverse set of challenges that represent
clinically occurring cases. Consequently, to develop a reliable method in
practice, we need not only to attain the desired accuracy on such a set, but
also to ensure that the method is capable of generalizing to novel sources of
anomalies.

By providing a downstream segmentation task, we remove any constraints on the
method design. One can use segmentation model features, its uncertainty
estimates, or an auxiliary model to detect outliers. Trained algorithms output
a single number called OOD score for every testing image; a higher score means
a higher outlier likelihood.

##  III Methods

###  III-A Methods selection

Several sub-topics, including anomaly detection (AD), novelty detection,
uncertainty estimation (UE), and outlier detection, share motivation and
methodology with OOD detection. Despite subtle differences between these
topics, the approaches are similar, and most of them can be applied to OOD
detection with minimal changes, as shown in [24]. So we follow the structure
of [24] and select the core methods from OOD detection, UE, and AD. In our
selection, we prioritize the methods already implemented for medical imaging,
e.g., in [11], [25], and [26].

As a universal baseline, the maximum probability of softmax outputs can be
used to detect OOD samples without any model modifications [4]. In practice,
however, the entropy of the same softmax outputs (Entropy) is used instead
[25, 11, 27]. We consider Entropy a starting point for all other approaches
and show its performance in our task.

The softmax entropy captures the total uncertainty, while the OOD measure
corresponds only to the epistemic uncertainty, as explained in [28]. Thereby,
one can use epistemic uncertainty estimation to improve over Entropy. Among
the others, Deep Ensemble [29] is considered the state-of-the-art approach to
UE. To use Ensemble, one computes mutual information or variance over several
predictions for a single image to obtain epistemic uncertainty map. An
alternative way to obtain multiple predictions is Monte-Carlo dropout (MCD)
[30], which we also include in our comparison.

Further, we include the approach of [11], which directly addresses OOD
detection on 3D medical images. The authors apply singular value decomposition
(SVD) to the network features and use singular values as an image embedding.
OOD score is calculated as the distance from a sample’s embedding to its
nearest neighbor from a training set.

A better uncertainty estimation can be obtained by modifying the downstream
model, although such modifications can harm the model’s performance. We
include one such popular modification, generalized ODIN (G-ODIN) [31], in our
study. Finally, OOD scores can be obtained with an auxiliary model, solely
dedicated to detecting anomalies in data. Such AD methods were extensively
compared in the Medical Out-of-Distribution (MOOD) challenge [12]. We choose
the best solution from MOOD 2022, implement, and include it in our experiments
under the name MOOD-1.

Discussing the auxiliary AD models, we intentionally exclude the
reconstruction-based methods (e.g., auto-encoders, generative-adversarial
nets) from our consideration. Firstly, these methods performed substantially
worse in MOOD 2022 than self-supervised learning-based ones (e.g., MOOD-1)
[26]. [32] also showed that this type of methods scores far behind self-
supervised learning. And [33] highlighted the severe limitations of auto-
encoders applied to OOD detection in a similar setup. Given this critique, we
do not include the reconstruction-based approaches in our experiments.

So, we consider the following methods: Entropy, Ensemble, MCD, SVD, G-ODIN,
and MOOD-1. Since some of them are designed for the downstream classification
task, we detail their adaptation to segmentation below.

###  III-B Methods implementation

To preserve a fair comparison, we add only trivial and unavoidable
modifications. We also test (in preliminary experiments) any additional
component or a critical hyperparameter of every method and select the best
performing setting.

####  III-B1 Entropy

Our downstream task is binary segmentation, where the sigmoid function is
applied to the network’s outputs. We note that two-classes softmax can be
derived from the sigmoid. Then, Entropy follows the implementation from [27]
and [11], computing the average entropy value over the predicted area (i.e.,
positive class). We set OOD score to 0 in the case of empty predicted mask.

####  III-B2 Ensemble

We train 5 U-Net models with different initializations and calculate the
uncertainty map as the voxel-wise standard deviation of the five corresponding
predictions. OOD score is the average of this uncertainty map.

####  III-B3 MCD

We implement MCD by introducing a dropout layer before every down- and up-
sampling layer in U-Net. To obtain an uncertainty map, we calculate voxel-wise
standard deviations of 5 inference steps with a dropout rate of 0.1. OOD score
is the average of this uncertainty map.

####  III-B4 SVD

We follow [11] without any changes.

####  III-B5 G-ODIN

We preserve the original structure of the G-ODIN output layer [31]; the only
difference is that we substitute the linear layers with the convolution ones.
These convolution layers have kernels of size $1^{3}$, so the procedure
remains equal to classification of every voxel. To obtain the uncertainty map,
we use the best reported G-ODIN DeConf-C∗ variant.

In Ensemble, MCD, and G-ODIN, using mutual information or entropy and
averaging the uncertainty map over only the predicted area, as in Entropy,
harms the performance. Thus, we use simple averaging of uncertainty.

####  III-B6 MOOD-1

The top-performing MOOD solutions generate synthetic anomalies and train a
network to segment them [26]. So our MOOD-1 implementation is based on this
cut-paste-segment approach, which won MOOD 2021 [34]. We then supplement it
with technical improvements from the 2022’s best solution, such as one-cycle
learning and ensembling over 5 models. The subject-level OOD score is
calculated as the mean of the top-100 anomaly probabilities.

####  III-B7 Volume predictor

To demonstrate that some semantic differences might be trivial from the
model’s perspective but not captured by other methods, we use the total volume
of prediction (positive class) as an OOD score. Since a predicted volume can
vary in any direction, we consider sample an outlier if the volume is below
$\frac{q}{2}$-th or above $100-\frac{q}{2}$-th percentile of the ID, thus
retaining $100-q$ TPR.

###  III-C Intensity Histogram Features

Figure 3: The proposed OOD detection method, called Intensity Histogram
Features (IHF). It consists of three steps: calculating a $m$-dimensional
vector as a histogram bin values from the preprocessed image (Step 1), fitting
and applying PCA to the occuring data, and calculating Mahalanobis distance
between a test vector and ID samples distribution (Step 3). We apply IHF to
the 3D images and illustrate the process using 2D axial slices for simplicity.
(*PCA is fitted once on all training data.)

To contest the DL algorithms, we propose an unsupervised method based on image
intensity histograms as embeddings. Our design is motivated by two other
works. Firstly, [11] showed that SVD can efficiently reduce full-image-sized
network features. We note a space for improvement in their method – one can
optimize the choice of the network’s layer to apply SVD. Here, [35] suggested
that the earlier network’s layers contain the most domain-specific
information. Following the latter suggestion, we hypothesize that we can
extract enough domain-specific information directly from the image (i.e., the
zeroth network’s layer). A histogram is a convenient way to do so.

We present our method, called Intensity Histogram Features (IHF),
schematically in Fig. 3. It consists of three steps: (1) calculating intensity
histograms of images and using them as vectors, (2) reducing their
dimensionality with PCA, and (3) running the outlier detection algorithms on
these vectors.

##### Step 1: preprocessing and histograms

All images undergo the same preprocessing pipeline to standardize the
intensity distribution:

  1. 1.

We interpolate images to the same spacing. So in all CT and MRI experiments,
we use $1\times 1\times 1.5$ mm.

  2. 2.

We clip image intensities to $[-1350;300]$ Hounsfield units for CT (a standard
lung window) and $[1^{th}$percentile; $99^{th}$percentile$]$ for MRI.

  3. 3.

We min-max-scale image intensities to the $[0,1]$ range.

Given a preprocessed image $x$, we then compute a probability density function
of its intensities in $m$ bins, a histogram $e(x)\in\mathbb{R}^{m}$, and
further use these vectors $e(x)$.

##### Step 2: Principal Component Analysis (PCA)

As an optional step, we use PCA to reduce the dimensions $m$. The main reason
to use it is that some outlier detection algorithms at Step 3 behave unstable
in high dimensional spaces. For instance, calculating Mahalanobis distance
require reversing the empirical sample covariance matrix, and this matrix is
likely to become ill-conditioned or singular with larger $m$.

Therefore, we fit PCAv once on the training data $E_{tr}$ to preserve
$v=99.99\%$ of the explained variance. This mostly eliminates the potential
instability while preserving the distribution properties. $E_{tr}$ consists of
row-vectors $e(x_{tr})$ for all training images $x_{tr}\in X_{tr}$. Further,
we use transformed vectors $\tilde{e}(x)=\text{PCA}_{v}(e(x))$.

##### Step 3: OOD detection algorithm

To calculate the OOD score for $x$, we can apply any distance- or density-
based outlier detection method. As in [36], we can calculate Mahalanobis
distance $S_{Mah}(x)$:

$S_{Mah}(x)=\sqrt{\left(\tilde{e}(x)-\hat{\mu}\right)^{T}\hat{\Sigma}^{-1}\left(\tilde{e}(x)-\hat{\mu}\right)},$
(1)

where $\hat{\mu}$ and $\hat{\Sigma}$ are the estimated mean and covariance
matrix on the training set, $\hat{\mu}=\frac{1}{|X_{tr}|}\sum_{x_{tr}\in
X_{tr}}\tilde{e}\left(x_{tr}\right)$ and
$\hat{\Sigma}=\frac{1}{|X_{tr}|}\sum_{x_{tr}\in
X_{tr}}\left(\tilde{e}(x_{tr})-\hat{\mu}\right)\left(\tilde{e}(x_{tr})-\hat{\mu}\right)^{T}$.
Alternatively, one can calculate the distance to the nearest neighbor (min-
distance) $S_{NN}(x)$, as in [11]:

$S_{NN}(x)=\min_{x_{tr}\in X_{tr}}||\tilde{e}(x)-\tilde{e}(x_{tr})||_{2}.$ (2)

Using $S_{Mah}$ (Eq. 1) and $S_{NN}$ (Eq. 2) corresponds to the methods called
IHF-Mah and IHF-NN, respectively. We include them in comparison and ablation
study independently.

##  IV Experimental setup

###  IV-A Downstream task

We have 3D CT and MRI datasets with a binary segmentation task. So we adhere
the standard approaches to train a segmentation model in all methods that
require the latter.

#### Data preprocessing

We describe preprocessing in IHF, Step 1 (Sec. III-C); it is the same in all
experiments and it is the minimum allowing the correct DL model training.

#### Architecture and training

In all experiments, we use 3D U-Net [37], a standard architecture for
segmentation. We train it on patches of $64$ axial slices, with a batch size
of $3$, Adam optimizer, and learning rate of $10^{-4}$ for $30$ epochs, $1000$
iterations each. In a batch, patches from different images are padded if
necessary. We minimize the sum of Binary Cross-Entropy and Focal Tversky
losses [38] to achieve high segmentation sensitivity.

##### Segmentation evaluation

We train all models on the training part of ID datasets. Then, we can evaluate
their segmentation quality on the corresponding testing part of the OOD
datasets, showing its possible decline. These segmentation results are given
in Tab. III for the CT and in Tab. IV for the MRI datasets.

###  IV-B OOD detection evaluation

Given the testing part of the ID dataset, we measure the OOD detection quality
against it for all the suggested OOD setups, similarly to the classification
task. Outliers occur rarely in practice, so we aim to measure detection
quality when most of the ID samples are being preserved w.r.t. relatively rare
OOD events. In this case, one of the most convenient classification metrics to
use is false positive rate at 95% true positive rate (FPR), so we consider FPR
our primary metric Nonetheless, for the consistency with other studies, we
report AUROC in the supplementary materials.

##  V Results

In this section, we report on our experiments and results. We start by
benchmarking all considered methods, then present the analysis of the
benchmark design and finish with the ablation study of the methods on
synthetic data.

TABLE I: Comparison of the considered OOD detection methods in terms of
FPR@TPR95% scores (lower is better). We highlight the best scores in every row
in bold and ranked the methods by their average performance. The first and
second sections correspond to CT and MRI setups, respectively.

OOD setup | IHF-NN | SVD | IHF-Mah | MOOD-1 | G-ODIN | Volume | MCD | Ensemble | Entropy
---|---|---|---|---|---|---|---|---|---
Location (Head) | .00 | .00 | .00 | .12 | .55 | .53 | .36 | .51 | .56
Location (Liver) | .51 | .13 | .64 | .56 | .56 | .84 | .89 | .93 | .78
Population (COVID-19) | .54 | .75 | .72 | .51 | .54 | .82 | .58 | .58 | .87
Scanner | .88 | .89 | .85 | .73 | .92 | .86 | .89 | .90 | .83
Synthetic (Elastic) | .15 | .37 | .67 | .16 | .59 | .81 | .42 | .37 | .84
Synthetic (Image noise) | .49 | .37 | .62 | .11 | .89 | .85 | .87 | .82 | .81
Population (Glioblastoma) | .00 | .00 | .00 | .10 | .21 | .01 | .85 | .81 | .86
Population (Healthy) | .00 | .00 | .00 | .11 | .00 | .00 | .88 | 1.0 | .85
Scanner | .00 | .00 | .00 | .15 | .00 | .74 | .63 | .66 | .89
Synthetic (K-space noise) | .00 | .36 | .00 | .88 | .88 | .90 | .82 | .77 | .73
Synthetic (Anisotropy) | .09 | .20 | .05 | .57 | .88 | .93 | .77 | .77 | .81
Synthetic (Motion) | .00 | .58 | .00 | .73 | .93 | .94 | .85 | .88 | .91
Synthetic (Image noise) | .47 | .33 | .47 | .30 | .56 | .71 | .78 | .75 | .75
CT average | .43 | .42 | .58 | .36 | .67 | .79 | .67 | .68 | .78
MRI average | .08 | .21 | .07 | .41 | .50 | .60 | .80 | .81 | .83

###  V-A Benchmarking

Tab. I presents the primary results of our study. Uncertainty-based methods,
the ones that are not designed for segmentation, mostly fail in the suggested
challenges. Entropy, Ensemble, MCD, and G-ODIN gives substantially worse
(higher) FPR than the other methods, with only G-ODIN slightly surpassing a
simple Volume predictor. Methods dedicated to segmentation perform better on
average. For instance, our implementation of the best MOOD 2022 solution,
MOOD-1, achieves $0.36$ and $0.41$ average FPR on CT and MRI data,
respectively. SVD improves further; it appears the only reliable studied
method, providing $0.42$ and $0.21$ mean FPR.

Then, SVD performance is contested by the proposed IHF. In a combination with
min-distance, IHF-NN provides the best average score across studied
challenges: 0.43 and 0.08 FPR, respectively. In a combination with Mahalanobis
distance, IHF-Mah provides practically worse results in the CT setups.
Although IHF-Mah performs weaker, it was historically the first and we
submitted it to MOOD 2022 ($m=150$, no PCA). We placed second (team
AIRI333http://medicalood.dkfz.de/web/) with the earliest IHF version,
supporting its robustness by the independent evaluation.

We also conduct an ablation study to verify IHF robustness. As shown in Fig.
4, we test IHF performance varying its only two parameters, the number of bins
($m$) and explained variance ratio ($v$). Our findings indicate a consistent
behavior and numerical stability regardless of the parameter choice, with a
slight trend of improved quality at a larger $m$.

Both IHF variants perform comparable or better on average than SVD and,
consequently, the other studied methods. Therefore, we conclude that the
histograms of image intensities are descriptive enough to detect most of the
suggested OOD cases, while neural networks might omit important domain
features that could be used in this problem. We thus hypothesize that neural
networks-based OOD detection can be further improved and leave this promising
direction for future research.

We present the same comparison in terms of AUROC in Tab. V. Although AUROC is
not our primary metric, it roughly preserves the same relative ranking of the
studied methods, not contradicting to our main message.

Figure 4: Dependence of IHF on its two hyperparameters: the number of
histogram bins ($m$) and explained variance in PCA ($v$). We give the results
for both IHF variants and CT and MRI setups.

###  V-B In-depth benchmark analysis

Further, we emphasize the significance of constructing a correct benchmark to
study the methods. Analysis of our experimental results suggests the
following:

  * •

The benchmark should contain diverse OOD challenges.

  * •

Challenges should represent clinically occurring cases.

  * •

Potential biases in the benchmark should be explored using simple methods,
such as IHF or Volume predictor.

It is often possible to develop a method tailored to specific OOD sources
where it thrives but fails in the other setups. For example, G-ODIN
demonstrates near-perfect results in the Population and Scanner setups on MRI
data but yields the worst scores in the others. In practice however, the
precise anomaly source is always unknown, and a general OOD detector, with an
acceptable average performance, is needed. The true method effectiveness can
be estimated only in the context of diverse challenges.

Secondly, OOD sources should accurately represent or simulate the clinically
occurring cases. For instance, the Synthetic (Noise) setup, as introduced in
[26] and reproduced in our study, is not supported by any medical imaging
process. MOOD-1 achieves the highest performance in this setup because its
training objective closely aligns with the anomaly synthesis process. However,
performing well in this and similar cases is of no clinical value and,
consequently, biases the methods’ evaluation towards explicitly unrealistic
scenarios.

Finally, our analysis reveals that OOD challenges might contain implicit but
trivial features. If a benchmark focuses solely on any such feature, we can
design a method that exploits this feature, leading to deceptive conclusions
about the generalized performance. Instead, we suggest using such methods to
reveal biased features beforehand. For example, near-perfect IHF results in
several setups demonstrate that certain anomalies are actually trivial
intensity changes, reinforcing the need of designing diverse benchmarks.

TABLE II: Fechner correlation coefficients between IHF-NN, Volume and the
other studied methods’ performance.

| IHF-NN | SVD | IHF-M | MOOD | GODIN | Vol | MCD | Ens | Ent
---|---|---|---|---|---|---|---|---|---
IHF-NN | 1.00 | 0.38 | 0.85 | 0.23 | -0.08 | 0.23 | 0.08 | -0.08 | -0.23
Volume | 0.23 | 0.54 | 0.38 | 0.38 | 0.38 | 1.00 | 0.23 | 0.08 | 0.23

To ensure methods’ generalization, we calculate the Fechner correlation
between their results and results of the simple methods. We show that, apart
from SVD, the other methods do not exhibit a strong correlation with the
Volume or IHF scores (Tab. II). So the examined methods mostly do not rely on
trivial features, such as image intensity distribution. However, SVD shows a
correlation of 0.54 with the Volume scores, suggesting its hindered
generalization on new sources of data with small difference in the predicted
area volume.

Figure 5: FPR under synthetically distorted data for every distortion severity
level. Blue line indicates method’s average trend across presented challenges
with 95% confidence interval. The other UE methods (MCD, Ensemble, and G-ODIN)
are excluded since their average trend is similar to Entropy.

###  V-C Ablation study on synthetic data

In Fig. 5, we show the OOD detection results on synthetic data for different
levels of distortion. A larger severity level indicates an image corruption
with the larger magnitude, where the magnitudes are chosen perceptually from 1
(barely noticeable distortion) to 5 (heavily distorted image). The general
trend is that more distorted images are easier to detect. Here, SVD exhibits
the steepest average slope and behaves almost linearly with the increasing
severity level, suggesting that we have considered challenging but solvable
tasks.

Different methods exhibit different sensitivity to the level of distortion
required for detection. Entropy and the other UE methods begin to operate
effectively only at level 3, while IHF can detect anomalies at the minimal
level. So we conclude that the methods should be studied across a wide range
of the anomaly severity levels. Additionally, we show that MOOD-1 depends more
on the OOD source than severity level: it completely fails in Motion and
K-space setups, while almost perfectly detects Noise and Elastic deformations
independently of the severity level. Moreover, MOOD-1 and IHF behave inversely
to each other in Noise and Elastic setups. Such diverse behavior suggests the
need to study methods also across a wide range of anomaly types.

##  VI Discussion

Besides benchmarking the OOD detection methods, our study also suggests
practical ideas on building a correct benchmark. We mainly highlight the data
diversity in multiple dimensions and clinical relevance of the setups.
However, we leave several critical questions undeveloped, thus opening the
corresponding research directions, which we discuss below.

####  VI-1 Benchmark design

We select the datasets that represent the clinically occurring sources of OOD
data. And we confirm the importance of the constructed challenges by the
degraded downstream performance (Tab. III, Tab. IV). However, we cannot state
with certainty that the highlighted sources are the only differences between
the distributions. While we name the primary difference in each case (e.g.,
acquisition parameters, patient populations), the distributions might differ
due to several causes and the other causes exist outside of our consideration.
So the development of a more refined benchmark with the controlled distinction
between the OOD sources would greatly facilitate further research in this
area.

Furthermore, the question of whether any semantic difference should be
considered abnormal requires additional investigation. For instance,
Population (Healthy) setup is considered OOD due to this apparent semantic
difference – healthy cases instead of pathology ones. Segmentation models
often yield correct empty predictions for such images regardless any OOD
detection decision. Rejecting a correct prediction in this case should be
considered a false negative instead of true positive outcome, lowering the
false positive rate.

####  VI-2 Uncertainty estimation

Our study verifies that the epistemic uncertainty is better suited for OOD
detection than the total uncertainty, since MCD and Ensemble achieve better
results than Entropy. Nevertheless, a question of how to aggregate the
uncertainty map into a single score remains open. On the one hand, aggregating
uncertainty over the predicted volume offers certain advantages compared to
the whole volume, especially when dealing with 3D images, where the area of
interest may occupy only a small portion of the entire image. While this
aggregation shows improved results for the Entropy method, it cannot rank
samples with an empty predicted mask and does not trigger on anomalies outside
of the predicted area. Contrary, aggregating uncertainty as a simple average
over the whole image provides better results for MCD and Ensemble. Developing
a reasonable UE method for 3D images, thereby, is a possible direction for
future research.

####  VI-3 Other IHF applications

Additionally, we explored alternative applications of the proposed IHF method.
We noted its strong performance in such medical imaging tasks as contrast
detection and domain classification. In this paper, however, we do not delve
into the potential IHF applications, as they lie beyond the scope of the OOD
detection problem.

##  VII Conclusion

In this paper, we have conducted an extensive investigation of OOD detection
on 3D medical images. Our results revealed that the established approaches,
including uncertainty estimation and anomaly detection ones, are unable to
provide reliable performance. They provide unacceptably high number of false
positives (0.31 mean FPR at best) and fail to generalize. They can also be
improved. To demonstrate this, we developed a histogram-based method, IHF,
that achieves comparable and often superior results to its competitors.
Thereby, we indicated that the distribution shifts in 3D medical imaging can
often be detected using intensity histograms, while the DL algorithms neglect
this domain feature. Although IHF achieves better average results, its
performance is surpassed in multiple challenges, emphasizing the need and
possibility for developing a robust and general OOD detection method.

We constructed and released the corresponding challenges as a benchmark for
OOD detection on 3D medical images, proposing IHF as a solid baseline to
contest new methods.

##### Acknowledgments

We acknowledge the National Cancer Institute and the Foundation for the
National Institutes of Health, and their critical role in the creation of the
free publicly available LIDC/IDRI Database used in this study.

## References

  * [1] M. Wang and W. Deng, “Deep visual domain adaptation: A survey,” _Neurocomputing_ , vol. 312, pp. 135–153, 2018.
  * [2] B. Kompa, J. Snoek, and A. L. Beam, “Second opinion needed: communicating uncertainty in medical machine learning,” _NPJ Digital Medicine_ , vol. 4, no. 1, p. 4, 2021.
  * [3] J. Yang, K. Zhou, Y. Li, and Z. Liu, “Generalized out-of-distribution detection: A survey,” _arXiv preprint arXiv:2110.11334_ , 2021.
  * [4] D. Hendrycks and K. Gimpel, “A baseline for detecting misclassified and out-of-distribution examples in neural networks,” _arXiv preprint arXiv:1610.02136_ , 2016.
  * [5] D. Hendrycks, S. Basart, M. Mazeika, M. Mostajabi, J. Steinhardt, and D. Song, “Scaling out-of-distribution detection for real-world settings,” _arXiv preprint arXiv:1911.11132_ , 2019.
  * [6] A. Mahmood, J. Oliva, and M. Styner, “Multiscale score matching for out-of-distribution detection,” _arXiv preprint arXiv:2010.13132_ , 2020\.
  * [7] A. G. Pacheco, C. S. Sastry, T. Trappenberg, S. Oore, and R. A. Krohling, “On out-of-distribution detection algorithms with deep neural skin cancer classifiers,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops_ , 2020, pp. 732–733.
  * [8] C. Berger, M. Paschali, B. Glocker, and K. Kamnitsas, “Confidence-based out-of-distribution detection: a comparative study and analysis,” in _Uncertainty for Safe Utilization of Machine Learning in Medical Imaging, and Perinatal Imaging, Placental and Preterm Image Analysis_. Springer, 2021, pp. 122–132.
  * [9] T. Cao, C.-W. Huang, D. Y.-T. Hui, and J. P. Cohen, “A benchmark of medical out of distribution detection,” _arXiv preprint arXiv:2007.04250_ , 2020\.
  * [10] G. Litjens, T. Kooi, B. E. Bejnordi, A. A. A. Setio, F. Ciompi, M. Ghafoorian, J. A. Van Der Laak, B. Van Ginneken, and C. I. Sánchez, “A survey on deep learning in medical image analysis,” _Medical image analysis_ , vol. 42, pp. 60–88, 2017.
  * [11] D. Karimi and A. Gholipour, “Improving calibration and out-of-distribution detection in deep models for medical image segmentation,” _IEEE Transactions on Artificial Intelligence_ , 2022.
  * [12] D. Zimmerer, J. Petersen, G. Köhler, P. Jäger, P. Full, K. Maier-Hein, T. Roß, T. Adler, A. Reinke, and L. Maier-Hein, “Medical Out-of-Distribution Analysis Challenge 2022,” Mar. 2022. [Online]. Available: https://doi.org/10.5281/zenodo.6362313
  * [13] B. Lambert, F. Forbes, S. Doyle, A. Tucholka, and M. Dojat, “Improving uncertainty-based out-of-distribution detection for medical image segmentation,” _arXiv preprint arXiv:2211.05421_ , 2022.
  * [14] S. G. Armato III, G. McLennan, L. Bidaut, M. F. McNitt-Gray, C. R. Meyer, A. P. Reeves, B. Zhao, D. R. Aberle, C. I. Henschke, E. A. Hoffman _et al._ , “The lung image database consortium (lidc) and image database resource initiative (idri): a completed reference database of lung nodules on ct scans,” _Medical physics_ , vol. 38, no. 2, pp. 915–931, 2011.
  * [15] S. Morozov, V. Gombolevskiy, A. Elizarov, M. Gusev, V. Novik, S. Prokudaylo, A. Bardin, E. Popov, N. Ledikhova, V. Chernina _et al._ , “A simplified cluster model and a tool adapted for collaborative labeling of lung cancer ct scans,” _Computer Methods and Programs in Biomedicine_ , vol. 206, p. 106111, 2021.
  * [16] E. B. Tsai, S. Simpson, M. P. Lungren, M. Hershman, L. Roshkovan, E. Colak, B. J. Erickson, G. Shih, A. Stein, J. Kalpathy-Cramer _et al._ , “The rsna international covid-19 open radiology database (ricord),” _Radiology_ , vol. 299, no. 1, pp. E204–E213, 2021.
  * [17] P. Bilic, P. F. Christ, E. Vorontsov, G. Chlebus, H. Chen, Q. Dou, C.-W. Fu, X. Han, P.-A. Heng, J. Hesser _et al._ , “The liver tumor segmentation benchmark (lits),” _arXiv preprint arXiv:1901.04056_ , 2019.
  * [18] M. Hssayeni, M. Croock, A. Salman, H. Al-khafaji, Z. Yahya, and B. Ghoraani, “Computed tomography images for intracranial hemorrhage detection and segmentation,” _Intracranial Hemorrhage Segmentation Using A Deep Convolutional Model. Data_ , vol. 5, no. 1, 2020.
  * [19] J. Shapey, A. Kujawa, R. Dorent, G. Wang, A. Dimitriadis, D. Grishchuk, I. Paddick, N. Kitchen, R. Bradford, S. R. Saeed _et al._ , “Segmentation of vestibular schwannoma from mri, an open annotated dataset and baseline algorithm,” _Scientific Data_ , vol. 8, no. 1, pp. 1–6, 2021\.
  * [20] R. Dorent, A. Kujawa, S. Cornelissen, P. Langenhuizen, J. Shapey, and T. Vercauteren, “Cross-Modality Domain Adaptation Challenge 2022 (crossMoDA),” May 2022. [Online]. Available: https://doi.org/10.5281/zenodo.6504722
  * [21] S. R. van der Voort, F. Incekara, M. M. Wijnenga, G. Kapsas, R. Gahrmann, J. W. Schouten, H. J. Dubbink, A. J. Vincent, M. J. van den Bent, P. J. French _et al._ , “The erasmus glioma database (egd): Structural mri scans, who 2016 subtypes, and segmentations of 774 patients with glioma,” _Data in brief_ , vol. 37, p. 107191, 2021.
  * [22] R. Souza, O. Lucena, J. Garrafa, D. Gobbi, M. Saluzzi, S. Appenzeller, L. Rittner, R. Frayne, and R. Lotufo, “An open, multi-vendor, multi-field-strength brain mr dataset and analysis of publicly available skull stripping methods agreement,” _NeuroImage_ , vol. 170, pp. 482–494, 2018.
  * [23] F. Pérez-García, R. Sparks, and S. Ourselin, “Torchio: a python library for efficient loading, preprocessing, augmentation and patch-based sampling of medical images in deep learning,” _Computer Methods and Programs in Biomedicine_ , p. 106236, 2021. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0169260721003102
  * [24] J. Yang, P. Wang, D. Zou, Z. Zhou, K. Ding, W. Peng, H. Wang, G. Chen, B. Li, Y. Sun _et al._ , “Openood: Benchmarking generalized out-of-distribution detection,” _arXiv preprint arXiv:2210.07242_ , 2022.
  * [25] A. Jungo and M. Reyes, “Assessing reliability and challenges of uncertainty estimations for medical image segmentation,” in _International Conference on Medical Image Computing and Computer-Assisted Intervention_. Springer, 2019, pp. 48–56.
  * [26] D. Zimmerer, P. M. Full, F. Isensee, P. Jäger, T. Adler, J. Petersen, G. Köhler, T. Ross, A. Reinke, A. Kascenas _et al._ , “Mood 2020: A public benchmark for out-of-distribution detection and localization on medical images,” _IEEE Transactions on Medical Imaging_ , vol. 41, no. 10, pp. 2728–2738, 2022.
  * [27] A. Mehrtash, W. M. Wells, C. M. Tempany, P. Abolmaesumi, and T. Kapur, “Confidence calibration and predictive uncertainty estimation for deep medical image segmentation,” _IEEE Transactions on Medical Imaging_ , vol. 39, no. 12, pp. 3868–3878, dec 2020. [Online]. Available: https://doi.org/10.1109%2Ftmi.2020.3006437
  * [28] L. Smith and Y. Gal, “Understanding measures of uncertainty for adversarial example detection,” 2018.
  * [29] B. Lakshminarayanan, A. Pritzel, and C. Blundell, “Simple and scalable predictive uncertainty estimation using deep ensembles,” _Advances in neural information processing systems_ , vol. 30, 2017.
  * [30] Y. Gal and Z. Ghahramani, “Dropout as a bayesian approximation: Representing model uncertainty in deep learning,” in _international conference on machine learning_. PMLR, 2016, pp. 1050–1059.
  * [31] Y.-C. Hsu, Y. Shen, H. Jin, and Z. Kira, “Generalized odin: Detecting out-of-distribution image without learning from out-of-distribution data,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2020, pp. 10 951–10 960.
  * [32] Y. Liang, J. Zhang, S. Zhao, R. Wu, Y. Liu, and S. Pan, “Omni-frequency channel-selection representations for unsupervised anomaly detection,” _arXiv preprint arXiv:2203.00259_ , 2022.
  * [33] F. Meissen, B. Wiestler, G. Kaissis, and D. Rueckert, “On the pitfalls of using the residual error as anomaly score,” _arXiv preprint arXiv:2202.03826_ , 2022.
  * [34] J. Cho, I. Kang, and J. Park, “Self-supervised 3d out-of-distribution detection via pseudoanomaly generation,” in _International Conference on Medical Image Computing and Computer-Assisted Intervention_. Springer, 2021, pp. 95–103.
  * [35] I. Zakazov, B. Shirokikh, A. Chernyavskiy, and M. Belyaev, “Anatomy of domain shift impact on u-net layers in mri segmentation,” in _International Conference on Medical Image Computing and Computer-Assisted Intervention_. Springer, Cham, 2021, pp. 211–220.
  * [36] K. Lee, K. Lee, H. Lee, and J. Shin, “A simple unified framework for detecting out-of-distribution samples and adversarial attacks,” _Advances in neural information processing systems_ , vol. 31, 2018.
  * [37] F. Isensee, P. Kickingereder, W. Wick, M. Bendszus, and K. H. Maier-Hein, “No new-net,” in _International MICCAI Brainlesion Workshop_. Springer, 2018, pp. 234–244.
  * [38] N. Abraham and N. M. Khan, “A novel focal tversky loss function with improved attention u-net for lesion segmentation,” in _2019 IEEE 16th international symposium on biomedical imaging (ISBI 2019)_. IEEE, 2019, pp. 683–687.

##  Appendix A Segmentation quality on ID and OOD datasets

TABLE III: Segmentation quality for all considered CT setups. We use Dice
score (DSC, higher is better) and average number of False Positives per image
(Avg. FP, lower is better) as our metrics.

Metric | ID dataset | OOD challenges
---|---|---
LIDC | Loc(Head) | Loc(Liver) | Scanner | Pop(COVID) | Syn(Noise) | Syn(Elastic)
DSC $(\uparrow)$ | $.50\pm.27$ | $.21\pm.41$ | $.03\pm.16$ | $.18\pm.24$ | $.00\pm.00$ | $.32\pm.30$ | $.13\pm.23$
Avg. FP $(\downarrow)$ | $5.2\pm 6.9$ | $2.6\pm 4.6$ | $10\pm 16$ | $6\pm 15$ | $13\pm 10$ | $12\pm 25$ | $5.4\pm 4.5$

TABLE IV: Segmentation quality for all considered MRI setups in terms of Dice
score (DSC).

Metric | ID dataset | OOD challenges |
---|---|---|---
VS-Seg | Pop(Healthy) | Pop(Glioma) | Scanner | Syn(K-space) | Syn(Anisotropy) | Syn(Motion) | Syn(Noise)
DSC $(\uparrow)$ | $.92\pm.05$ | $.85\pm.36$ | $.86\pm.34$ | $.89\pm.12$ | $.83\pm.023$ | $.87\pm.01$ | $.91\pm.06$ | $.57\pm.41$

##  Appendix B OOD detection results in terms of AUROC

TABLE V: Comparison of the considered OOD detection methods in terms of AUROC
scores (higher is better). We highlight the best scores in every row in bold
and ranked the methods by their average performance. The first and second
sections correspond to CT and MRI setups, respectively.

OOD setup | IHF-NN | IHF-Mah | SVD | G-ODIN | Volume | MOOD-1 | MCD | Ensemble | Entropy
---|---|---|---|---|---|---|---|---|---
Location (Head) | 1.0 | 1.0 | 1.0 | .83 | .73 | .83 | .85 | .79 | .62
Location (Liver) | .89 | .85 | .97 | .88 | .65 | .61 | .42 | .45 | .67
Population (COVID-19) | .88 | .83 | .74 | .86 | .76 | .66 | .79 | .80 | .72
Scanner | .73 | .73 | .58 | .72 | .68 | .51 | .58 | .55 | .65
Synthetic (Elastic) | .97 | .83 | .86 | .85 | .77 | .78 | .84 | .85 | .65
Synthetic (Image noise) | .81 | .77 | .84 | .75 | .67 | .80 | .56 | .61 | .59
Population (Glioblastoma) | 1.0 | 1.0 | 1.0 | .96 | .68 | .87 | .44 | .41 | .14
Population (Healthy) | 1.0 | 1.0 | 1.0 | 1.0 | .68 | .86 | .44 | .16 | .15
Scanner | 1.0 | 1.0 | 1.0 | 1.0 | .77 | .83 | .70 | .74 | .59
Synthetic (K-space noise) | 1.0 | 1.0 | .86 | .81 | .66 | .24 | .56 | .63 | .66
Synthetic (Anisotropy) | .98 | .98 | .94 | .81 | .68 | .57 | .63 | .63 | .71
Synthetic (Motion) | .99 | 1.0 | .75 | .78 | .68 | .48 | .57 | .54 | .57
Synthetic (Image noise) | .81 | .83 | .85 | .88 | .66 | .77 | .58 | .58 | .56
CT average | .88 | .84 | .83 | .82 | .71 | .70 | .67 | .68 | .65
MRI average | .97 | .97 | .92 | .89 | .69 | .66 | .56 | .53 | .48


$\lim_{n\rightarrow\infty}\mathbb{P}\left(\left|\left|D_{n}^{\frac{1}{2}}S_{ZZ}D_{n}^{\frac{1}{2}}-\Sigma_{Z,n}\right|\right|>\epsilon\right)=0$

and

(31)
$\lim_{n\rightarrow\infty}d_{BL}\left(\sqrt{n}\Sigma_{Z,n}^{-\frac{1}{2}}D_{n}^{\frac{1}{2}}S_{Z\epsilon}\Sigma^{-\frac{1}{2}},V\right)=0$

where $\textnormal{vec}(V)\sim\mathcal{N}(0,I)$.

###### Proof.

Throughout the proof we write matrices as $3\times 3$ block matrices such that
the top-left block is $r_{1}\times r_{1}$, the middle block is
$(r_{2}-r_{1})\times(r_{2}-r_{1})$, and the bottom-left block is
$(d-r_{2})\times(d-r_{2})$. We use a superscript to denote the block index,
e.g., $S_{ZZ}^{13}$ denotes the top-right block of $S_{ZZ}$. Let
$\tilde{Z}_{t}=\sum_{s=1}^{t}(1-n^{-\beta})^{t-s}\epsilon_{s}$ and
$\psi_{t}=\sum_{s=1}^{t}(1-n^{-\beta})^{t-s}X_{s-1}$ so that

(32) $Z_{t}=\tilde{Z}_{t}+\left(\Gamma_{n}-I\right)\psi_{t}.$

Let $\Lambda_{n}$ be the diagonal matrix such that
$\Lambda_{n,ii}=1-|\Gamma_{n,ii}|$ for $1\leq i\leq d$. Then, since
$\Gamma_{n}$ and $I-\Lambda_{n}$ commute, with $C$ denoting some generic
constant not depending on $t$ or $n$,

(33)
$\begin{split}\mathbb{E}\left|\left|(I-\Lambda_{n})\psi_{t}\right|\right|^{2}&=\sum_{i,j=1}^{t}(1-n^{-\beta})^{2t-i-j}\text{tr}\left((I-\Lambda_{n})\mathbb{E}\left(X_{i-1}X_{j-1}^{T}\right)(I-\Lambda_{n})\right)\\\
&\leq C\sum_{i,j=1}^{t}(1-n^{-\beta})^{2t-i-j}\left|\left|\Gamma^{i\vee
j-i\wedge j}(I-\Lambda_{n})\right|\right|\\\ &\leq
C\sum_{i=0}^{t-1}(1-n^{-\beta})^{i}\sum_{j=0}^{i-1}\left|\left|\left((1-n^{-\beta})\Gamma\right)^{j}(I-\Lambda_{n})\right|\right|\end{split}$

where $a\vee b=\max\\{a,b\\}$ and $a\wedge b=\min\\{a,b\\}$ and the first
inequality follows from Lemma F.3 and the Cauchy-Schwartz inequality. This
inequality yields a result equivalent to equation (40) in [21]. In particular,
we deduce that

$\sup_{1\leq t\leq
n}\mathbb{E}\left|\left|(\Gamma_{n}-I)\psi_{t}\right|\right|^{2}=o(n)$

from which it follows that

(34) $S_{Z\epsilon}=S_{\tilde{Z}\epsilon}+o_{p}(1).$

We first prove (30). For ease of notation, we write

$S_{n}=D_{n}^{\frac{1}{2}}S_{ZZ}D_{n}^{\frac{1}{2}}.$

Since $D_{n}^{11}=n^{-\beta}I_{r_{1}}$, essentially the same proof as that of
Lemma 3.1.(iii) in [21] using (33) shows that

$S_{n}^{11}=n^{-\beta}S^{11}_{\tilde{Z}\tilde{Z}}+o_{p}(1)=\frac{1}{2}\Sigma^{11}+o_{p}(1),$

where the latter equality follows from Theorem 3.2 and the fact that
$n^{-\beta}\mathbb{E}(S_{\tilde{Z}\tilde{Z}})\rightarrow\Sigma/2$ for
$n\rightarrow\infty$. Similarly, the proof of Lemma 3.5.(ii) in [21] can be
adapted to show that

$S_{n}^{33}=\left(D_{n}^{33}\right)^{\frac{1}{2}}S^{33}_{XX}\left(D_{n}^{33}\right)^{\frac{1}{2}}+o_{p}(1)=\Sigma^{33}_{X,n}+o_{p}(1),$

where the latter equality follows from Theorem 3.2 and Lemma F.3 and
$\Sigma^{33}_{X,n}$ is defined as in Lemma F.3 but emphasizing the dependence
on $n$. For the middle block, using the recursive relations
$Z_{t}=(1-n^{-\beta})Z_{t-1}+\Delta X_{t}$ and $\Delta
X_{t}=(\Gamma_{n}-I)X_{t-1}+\epsilon_{t}$, we can write

$\left(1-\left(1-n^{-\beta}\right)^{2}\right)S^{22}_{ZZ}=S^{22}_{\Delta
XZ}+S^{22}_{Z\Delta X}+S^{22}_{\Delta X\Delta X}+o_{p}(1).$

It follows from Theorem 3.2, that $S^{22}_{\Delta X\Delta
X}=\Sigma^{22}+o_{p}(1)$. For the other two terms, we use (34) and write
$S^{22}_{Z\Delta
X}=S^{22}_{\tilde{Z}\epsilon}+S^{22}_{ZX}(\Gamma_{n}^{22}-I)^{T}+o_{p}(1)$.
Using the recursive relations and (34) once more yields

$\left(I-\left(1-n^{-\beta}\right)\Gamma^{22}_{n}\right)S^{22}_{XZ}=S^{22}_{X\epsilon}+S^{22}_{\epsilon\tilde{Z}}+S^{22}_{\epsilon\epsilon}+S^{22}_{XX}(\Gamma^{22}_{n}-I)+o_{p}(1)$

It follows from Theorem 3.2 that the first two terms tend to 0 in probability
for $n\rightarrow\infty$ and
$\Sigma^{22}_{\epsilon\epsilon}=\Sigma^{22}+o_{p}(1)$. By Assumption U.5, for
$n$ large enough and $i\leq r_{2}$,

$1-\left|\Gamma_{n,ii}\right|=\left|1-\left|\Gamma_{n,ii}\right|\right|\leq\left|1-\Gamma_{n,ii}\right|\leq
1-\left|\Gamma_{n,ii}\right|,$

i.e., $|1-\Gamma_{n,ii}|=1-|\Gamma_{n,ii}|$ and, in particular,
$\Gamma_{n,ii}\in\mathbb{R}$ and $\kappa\in\mathbb{R}$. Thus, if we define
$K\in\mathbb{R}^{(r_{2}-r_{1})\times(r_{2}-r_{1})}$ diagonal with
$K_{ii}=\kappa_{i}$ for $i\leq r_{2}-r$ and $K_{ii}=\kappa_{i}^{-1}$
otherwise, we get

$\displaystyle n^{\beta}\left(I-\Gamma^{22}_{n}\right)$
$\displaystyle\rightarrow K$ $\displaystyle
n^{\beta}\left(\Lambda_{n}^{22}-I\right)$ $\displaystyle\rightarrow K$
$\displaystyle n^{\beta}\left(I-(1-n^{-\beta})\Gamma_{n}^{22}\right)$
$\displaystyle\rightarrow K+I$

for $n\rightarrow\infty$. Lemma F.3 then yields

$\left(\Gamma_{n}^{22}-I\right)S_{XZ}^{22}=\left(K+I\right)^{-1}\left(K\Sigma^{22}+K^{\frac{1}{2}}\Sigma_{X,n}^{22}K^{\frac{1}{2}}\right)+o_{p}(1).$

and
$(I-\Lambda^{22}_{n})^{\frac{1}{2}}\otimes(I-\Lambda^{22}_{n})^{\frac{1}{2}}(I-\Gamma_{n}^{22}\otimes\Gamma_{n}^{22})\rightarrow(K\otimes
K)^{\frac{1}{2}}(I\otimes K+K\otimes I)^{-1}$ for $n\rightarrow\infty$ so that
(keeping in mind that $K_{ii}>0$ for all $i$)

$\Sigma^{22}_{X,n}\rightarrow
K^{\frac{1}{2}}\int_{0}^{\infty}e^{-sK}\Sigma^{22}e^{-sK}dsK^{\frac{1}{2}}=K^{\frac{1}{2}}\Omega^{22}K^{\frac{1}{2}}.$

Then, using the relation $\Sigma^{22}-\Omega^{22}K=K\Omega^{22}$, the limiting
expression simplifies to

$(\Gamma_{n}^{22}-I)S_{XZ}^{22}=(K+I)^{-1}K^{2}\Sigma^{22}.$

Finally, since $n^{\beta}(1-(1-n^{-\beta})^{2})=2+o(1)$ and
$D_{n}^{22}=n^{-\beta}I_{r_{2}-r_{1}}$, we find

$\displaystyle S_{n}^{22}$
$\displaystyle=\frac{1}{2}\left(\Sigma^{22}+\left(K+I\right)^{-1}K^{2}\Omega^{22}+\Omega^{22}K^{2}\left(K+I\right)^{-1}\right)+o_{p}(1)$
$\displaystyle=\frac{1}{2}\left((K+I)^{-1}K\Omega^{22}+\Omega^{22}K(K+I)^{-1}\right)+o_{p}(1)$
$\displaystyle=\frac{1}{2}(I+K)^{-1}\left(2K\Omega^{22}K+\Sigma\right)(I+K)^{-1}+o_{p}(1).$

We have yet to characterize the asymptotic behaviour of the off-diagonal
blocks. First, note that by (23) in [21] and (32), we get

$S_{ZZ}^{32}-S_{XZ}^{32}=-n^{-\beta}\left(S_{\psi\psi}^{32}(\Gamma_{n}^{22}-I)-S_{\psi\tilde{Z}}^{32}\right)$

so that (33) and Theorem 3.2 yield

$S_{n}^{32}-\left(D_{n}^{33}\right)^{\frac{1}{2}}S_{XZ}^{32}\left(D_{n}^{22}\right)^{\frac{1}{2}}=o_{p}(1).$

Similar to above, we have

$\left(I-(1-n^{-\beta})\Gamma_{n}^{33}\right)S_{XZ}^{32}=S^{32}_{XX}\left(\Gamma_{n}^{22}-I\right)+o_{p}(1).$

But then, because $n^{-\beta}(I-\Gamma_{n}^{33})=o(1)$, we find that

$n^{-\frac{\beta}{2}}\left(I-\Lambda_{n}^{33}\right)^{\frac{1}{2}}\left(I-\left(1-n^{-\beta}\right)\Gamma_{n}^{33}\right)^{-1}=o(1)$

and, by Lemma F.3 and Theorem 3.2, $S_{XX}^{32}(\Gamma_{n}^{22}-I)=o_{p}(1)$.
Thus, $S_{n}^{32}=o_{p}(1)$. A similar argument show that
$S_{n}^{31}=o_{p}(1)$ so that the only block left is $S_{n}^{21}$. As a
consequence of (32), we have

$n^{-\beta}\left|\left|S_{ZZ}^{12}-S_{\tilde{Z}Z}^{12}\right|\right|=n^{-\beta}\left|\left|(\Gamma_{n}^{11}-I)S_{\psi\psi}^{12}(\Gamma_{n}^{22}-I)+(\Gamma_{n}^{11}-I)S_{\psi\tilde{Z}}^{12}\right|\right|$

and arguments like the one employed in the proof of Lemma 3.1 in [21] in
combination with (33) show that the right hand side is $o_{p}(1)$. Using the
recursive relations for $Z_{t}$, $\tilde{Z}_{t}$, and $\Delta X_{t}$ in
combination with (34) and Theorem 3.2, we have

$\displaystyle\left(1-(1-n^{-\beta})^{2}\right)S_{\tilde{Z}Z}^{12}$
$\displaystyle=S^{12}_{\epsilon Z}+S_{\tilde{Z}\Delta
X}^{12}+S_{\epsilon\Delta X}^{12}+o_{p}(1)$
$\displaystyle=S_{\tilde{Z}X}^{12}(\Gamma_{n}^{22}-I)+S^{12}_{\epsilon\epsilon}+o_{p}(1).$

An application of Lemma F.3 and Theorem 3.2 yields

$S_{\tilde{Z}X}^{12}=\Omega^{12}_{n}K^{\frac{1}{2}}+o_{p}(1)$

where

$\Omega^{12}_{n}=\Sigma^{12}\left(I-(1-n^{-\beta})\Gamma_{n}^{22}\right)n^{-\frac{\beta}{2}}\left(I-\Lambda_{n}^{22}\right)^{-\frac{1}{2}}\rightarrow\Sigma^{12}K^{\frac{1}{2}}\left(I+K\right)^{-1}$

for $n\rightarrow\infty$. In conclusion, since
$D_{n}^{22}=n^{-\beta}I_{r_{2}-r_{1}}$ and $D_{n}^{33}=n^{-\beta}I_{r_{1}}$,
we get

$S_{n}^{12}=\frac{1}{2}\left(\Sigma^{12}+\Sigma^{12}K(I+K)^{-1}\right)+o_{p}(1)=\frac{1}{2}\Sigma^{12}+o_{p}(1).$

Collecting all the limiting expressions, we define

$\bar{K}=\begin{pmatrix}I_{r_{1}}&0\\\
0&K\end{pmatrix},\quad\Omega=\int_{0}^{\infty}e^{-s\bar{K}}\begin{pmatrix}\Sigma^{11}&\Sigma^{12}\\\
\Sigma^{21}&\Sigma^{22}\end{pmatrix}e^{-s\bar{K}}ds$

and observe that

$\displaystyle\Sigma^{11}$
$\displaystyle=\left((I+\bar{K})^{-1}\left(2\bar{K}\Omega\bar{K}+\Sigma\right)(I+\bar{K})^{-1}\right)^{11}$
$\displaystyle\Sigma^{12}$
$\displaystyle=\left((I+\bar{K})^{-1}\left(2\bar{K}\Omega\bar{K}+\Sigma\right)(I+\bar{K})^{-1}\right)^{12}.$

so that $S_{n}=\Sigma_{Z,n}+o_{p}(1)$ with

$\Sigma_{Z,n}=\frac{1}{2}\begin{pmatrix}(I+\bar{K})^{-1}\left(2\bar{K}\Omega\bar{K}+\Sigma\right)(I+\bar{K})^{-1}&0\\\
0&\Sigma^{33}_{X,n}\end{pmatrix}.$

To see that $\Sigma_{Z,n}$ is asymptotically invertible and bounded simply
note that $\bar{K}_{ii}\in[0,1]$ for all $1\leq i\leq r_{2}$ and $\Sigma$ is
positive definite. Therefore, the top left block of $\Sigma_{Z,n}$ is some
fixed positive definite matrix for all $n\in\mathbb{N}$. The result then
follows from Lemma F.3.

Once (30) has been established, the proof of (31) is completely analogous to
the proof of (9). ∎

###### Proof of Theorem 4.2.

Fix some $\beta\in(\frac{1}{2},1)$ and a sequence
$(\theta_{n})_{n\in\mathbb{N}}\subset\Theta$ so that
$\theta_{0,n}=(\Gamma_{0,n},\Sigma_{0,n},\cdot)=(F_{\theta_{n}}\Gamma_{n}F_{\theta_{n}}^{-1},F_{\theta_{n}}\Sigma_{n}F_{\theta_{n}}^{-1},\cdot)$
satisfies U. Now, by Lemma D.1 and D.2, there exists a sub-sequence
$(\theta_{n_{k}})_{k\in\mathbb{N}}$ such that

$\hat{t}^{2}_{IV,\Gamma_{n_{k}}}=nS_{\epsilon
Z}S_{ZZ}^{-1}S_{Z\epsilon}=nS_{\epsilon
Z}F_{\theta_{n_{k}}}^{T}\left(F_{\theta_{n_{k}}}S_{ZZ}F_{\theta_{n_{k}}}^{T}\right)^{-1}F_{\theta_{n_{k}}}S_{Z\epsilon}\rightarrow_{w}\chi^{2}_{d^{2}}$

for $k\rightarrow\infty$. This proves that
$\hat{t}^{2}_{IV,\Gamma_{n}}\rightarrow_{w}\chi^{2}_{d^{2}}$ for any sequence
$(\theta_{n})_{n\in\mathbb{N}}\subset\Theta$ from which uniform convergence
follows by Lemma 1 in [14]. To prove uniform asymptotic level of the
corresponding confidence region, note that the limiting distribution is
continuous with respect to Lebesgue measure so that a similar argument as
before works in this case. ∎

##  Appendix E Simulations

###  E.1. Confidence intervals

For each $n\in\\{50,75,100\\}$ and $d\in\\{3,4,5\\}$ the simulation experiment
is repeated 1000 times. In each repetition, for $i,j=1,...,d$, we draw
$U_{ij}\sim\text{Unif}([0,1])$ and set $\Gamma=U^{-1}\Lambda_{n}U$ where
$\Lambda_{n}\in\mathbb{R}^{d\times d}$ is diagonal with $\Lambda_{n,11}=1$ and
$\Lambda_{n,ii}=1-(1/n)^{1/(i-1)}$ for $i=2,...,d$. We then sample
$\epsilon_{t}\sim\mathcal{N}(0,\Sigma)$ i.i.d. for $t=1,...,n$ with

$\Sigma=\frac{1}{2}\left(I+\mathbbm{1}\mathbbm{1}^{T}\right)$

where $\mathbbm{1}=(1,...,1)^{T}\in\mathbb{R}^{d}$ and let

$X_{t}=\Gamma X_{t-1}+\epsilon_{t}\text{ for } t=1,...,n,\quad X_{0}=0.$

$X_{t}$ is a sample from a VAR(1) process with $\theta=(\Gamma,\Sigma,\cdot)$.
We then compute $CI_{b}$, $CI_{IV}$, and $CI_{LA}$ for this sample and record
the length of each confidence interval and whether it contains $\Gamma_{11}$.

###  E.2. Predictive regression testing

For both simulation experiments we fix $d=4$ and $\alpha=0.1$. This implies
that $\tilde{\Gamma}\in\mathbb{R}^{3\times 3}$. The two regimes correspond to
two different choices of $\tilde{\Gamma}$:

  * •

_Mixed Regime_ : In this setting $\tilde{\Gamma}$ is chosen as above, that is,
with roots of differing proximity to unity and with random eigenvectors
sampled anew for every simulation run.

  * •

_Non-stationary Regime_ : In this setting we set $\tilde{\Gamma}=I$ so that
$\tilde{X}_{t}$ is a random walk.

In both regimes the errors are i.i.d. Gaussian with covariance matrix $\Sigma$
as given above.

To obtain Figure 1, we do the following: For each $n\in\\{10,20,...,200\\}$,
we draw two samples $X_{t}$ from the VAR(1) processes given by the two choices
of $\tilde{\Gamma}$ and under the null $H_{0}:\gamma=0$. We then compute the
three tests on both samples recording whether the null was rejected or not.
This is repeated 1000 times and the rejection rate is the proportion of times
the null was rejected across all simulations.

For Figure 2, we do essentially the same thing except that we now fix $n=100$
and perform the experiment across different choices of $\gamma\neq 0$. In
particular, we run the experiment for $\gamma=\delta\mathbbm{1}$,
$\delta\in\\{0.005,0.01,...,0.1\\}$ and record the proportion of times the
null was rejected across all 1000 simulations.

###  E.3. EAM

We present a slightly generalized version of the EAM-algorithm. EAM stands for
Evaluation-Approximation-Maximization and the algorithm can more or less be
split into three steps. It is an algorithm for solving problems of the
following form

$\displaystyle\sup f(x)$ $\displaystyle\text{s.t. }g(x)\leq c(x)$

over $x\in\mathcal{X}\subset\mathbb{R}^{p}$ where $f$, $g$, and $c$ are fixed
scalar functions sufficiently smooth and satisfying certain requirements. In
[12] it is required that $f(x)=v^{T}x$ for some $v\in\mathbb{R}^{p}$. Here we
only require that $f(x)$ is convex and twice continuously differentiable on
$\mathcal{X}$. We assume that $c$ is costly to evaluate. Without going into
too much detail, the algorithm proceeds as follows:

  1. (1)

_Initialization:_ Randomly sample initial points $x^{(1)},...,x^{(k)}$ from
$\mathcal{X}$ and evaluate $c(x^{(i)})$ for $i=1,...,k$. Set $L=k$.

  2. (2)

Iterate the following three steps until convergence:

    1. (a)

_E-step_ : Evaluate $c(x^{(L)})$ and pick current optimum

$y^{*,L}=\max\left\\{f\left(x^{(i)}\right):g\left(x^{(i)}\right)\leq
c\left(x^{(i)}\right),i=1,...,L\right\\}$

    2. (b)

_A-step_ : Approximate $x\mapsto c(x)$ by a Gaussian process regression model,
with mean $\mu$, constant variance $\sigma^{2}$, and covariance kernel
$K_{\beta}(x-x^{\prime})=\exp(-\sum_{i=1}^{p}|x_{i}-x^{\prime}_{i}|^{2}/\beta_{i})$,
fitted on $(c(x^{(i)}),x^{(i)})$, $i=1,...,L$. Fitting the model yields the
mean function $c_{L}(x)$ and the variance function $s_{L}(x)$ as well as the
fitted parameters $\hat{\mu}_{L},\hat{\sigma}_{L},\hat{\beta}_{L}$.

    3. (c)

_M-step_ : With probability $1-\epsilon$, let

$x^{(L+1)}=\operatorname*{arg\,max}_{x\in\mathcal{X}}EI_{L}(x)$

and with probability $\epsilon$ draw $x^{(L+1)}$ randomly from $\mathcal{X}$.
Set $L=L+1$.

$EI_{L}(x)$ is the _expected improvement function_ and it is given by

$EI_{L}(x)=\left(f(x)-y^{*,L}\right)_{+}\left(1-\Phi\left(\frac{g(x)-c_{L}(x)}{\hat{\sigma}_{L}s_{L}(x)}\right)\right)$

where $\cdot_{+}=\max\\{\cdot,0\\}$ and $\Phi$ is the standard normal CDF.
Note that the optimization problem in the M-step can be reformulated as a
constrained optimization problem with smooth objective function and smooth
constraints for which all derivatives are known and it can therefore be solved
with standard solvers. A key observation is that we only evaluate $c$ once per
iteration. In practice this results in far fewer evaluations of $c$ when
compared to, say, grid methods.

This algorithm was intended to compute confidence intervals that arise as
projections of confidence regions exactly as is the case for $CI_{b}$. In this
case we would simply take $x=\Gamma$ and let $f(\Gamma)=\Gamma_{11}$,
$g(\Gamma)=\hat{t}^{2}_{\Gamma}$, and
$c(\Gamma)=\tilde{q}_{n,\Gamma}(\alpha)$. It does not really matter that
$\Gamma$ is a matrix since we can just vectorize it and redefine all the
functions correspondingly. This would give us the upper bound of $CI_{b}$ and
the lower bound can be found by taking $f(\Gamma)=-\Gamma_{11}$.

Similarly, the EAM-algorithm can be used to compute $\phi_{b}$. This is done
by letting $x=\tilde{\Gamma}$ and then
$f(\tilde{\Gamma})=\hat{t}^{2}_{0|\tilde{\Gamma}}$ with $g$ and $c$ as before,
but for $\tilde{\Gamma}$ instead of $\Gamma$. Note that in both cases $f$ and
$g$ are polynomials of $\text{vec}(x)$ and are therefore smooth with known
derivatives of all orders.222All code used for the simulations can be found at
https://github.com/cholberg/unif_inf_var. Our implementation of the EAM-
algorithm is based on [13].

##  Appendix F Supplementary results

Here we give some supplementary results that are useful throughout. We assume
for the entirity of this section that Assumptions M and U are true with the
restriction $F_{\theta}=I$. The first lemma collects some convergence results
related to the normalizing matrix, $H$.

######  Lemma F.1.

We have

  1. (a)

$\liminf_{n\rightarrow\infty}\inf_{\theta\in\Theta}\sigma_{min}(H)>0$.

  2. (b)

$\sup_{R_{n,d}}\sigma_{min}(H)^{-1}=O(n^{-\eta})$.

  3. (c)

$\sup_{R_{n,0}}\sup_{1\leq k\leq j\leq d}(1-|\lambda_{i_{j}}|)|H_{kj}|=O(1)$.

  4. (d)

$\sup_{R_{n,0}}||\Gamma^{t}||\leq C(1-\log(n)/n)^{t}$ for a constant $C\geq 0$
not depending on $\theta$ or $n$.

  5. (e)

$\sup_{R_{n,0}}||\sum_{t=0}^{n-2}\Gamma^{t}\Sigma(\Gamma^{t})^{T}||=O(n/\log
n)$

  6. (f)

$\sup_{R_{n,0}}||\sum_{t=0}^{n-2}t\Gamma^{t}\Sigma(\Gamma^{t})^{T}||=O(n/\log
n)$

###### Proof.

We start with the proof of (a). We have

(35)
$H=\frac{1}{n}\sum_{t=1}^{n}\sum_{s=1}^{t-1}\Gamma^{t-1-s}\Sigma\left(\Gamma^{t-1-s}\right)^{T}.$

Since every matrix in the summand is positive semidefinite, we get

$\inf_{\theta\in\Theta}\sigma_{min}(H)\geq\inf_{\theta\in\Theta}\frac{1}{n}\sum_{t=1}^{n-1}\sigma_{min}(\Sigma)=\frac{n-1}{n}\inf_{\theta\in\Theta}\sigma_{min}(\Sigma)>\frac{n-1}{n}c$

where $c=\inf_{\theta\in\Theta}\sigma_{min}(\Sigma)>0$ by Assumption U.2.

For the proof of (b), note that, for any $n$ large enough and $\theta\in
R_{n,d}$, $\Gamma$ is diagonal and therefore we have from (35)

$\displaystyle\sigma_{min}(H)$
$\displaystyle\geq\frac{\sigma_{min}(\Sigma)}{n}\sum_{t=1}^{n-1}\sum_{s=0}^{t-1}|\lambda_{min}(\Gamma)|^{2s}$
$\displaystyle\geq\frac{\sigma_{min}(\Sigma)}{n}\sum_{t=1}^{n-1}\frac{1-(1-n^{-\eta})^{2t}}{1-(1-n^{-\eta})^{2}}$
$\displaystyle=\frac{\sigma_{min}(\Sigma)}{n(1-(1-n^{-\eta})^{2})}\left(n-1-\sum_{t=1}^{n-1}(1-n^{-\eta})^{2t}\right)$
$\displaystyle=\frac{\sigma_{min}(\Sigma)}{n(1-(1-n^{-\eta})^{2})}\left(n-\frac{1-(1-n^{-\eta})^{2n}}{1-(1-n^{-\eta})^{2}}\right).$

For $n_{0}$ large enough, we get

$\frac{1}{n}\left(n-\frac{1-(1-n^{-\eta})^{2n}}{1-(1-n^{-\eta})^{2}}\right)\geq\frac{1}{2},\quad\forall
n\geq n_{0}$

and, thus, with $C=2\sup_{\theta\in\Theta}\sigma_{min}(\Sigma)^{-1}<\infty$,
we have that, for all $n\geq n_{0}$,

$\sup_{\theta\in R_{n,d}}\sigma_{min}(H)^{-1}\leq
C\left(1-\left(1-n^{-\eta}\right)^{2}\right).$

Since $n^{\eta}(1-(1-n^{-\eta})^{2})\rightarrow 2$ for $n\rightarrow\infty$
and all $\eta>0$, this proves (b).

Turning to the proof of (c), we immediately see that the statement is true
when $|\lambda_{1}| \leq 1-\alpha$ so we can assume without loss of generality
that all eigenvalues are strictly greater than $1-\alpha$ and, thus, that
$\Gamma$ is diagonal. We then get, for $1\leq k\leq j\leq d$,

$|H_{kj}|\leq\frac{|\Sigma_{kj}|}{n}\sum_{t=1}^{n-1}\sum_{s=0}^{t-1}|\lambda_{i_{k}}\lambda_{i_{j}}|^{s}\leq\frac{|\Sigma_{kj}|}{n}\sum_{t=1}^{n-1}\sum_{s=0}^{t-1}|\lambda_{i_{j}}|^{s}=\frac{|\Sigma_{kj}|}{n(1-|\lambda_{i_{j}}|)}\left(n-\frac{1-|\lambda_{i_{j}}|^{n}}{1-|\lambda_{i_{j}}|}\right)$

where the second inequality follows from the fact that $|\lambda_{i}|\leq 1$
for all $1\leq i\leq N$. But then (c) follows from the fact that
$1-|\lambda_{1}|\geq\log(n)/n$ for any $\theta\in R_{n,0}$.

We now prove (d). For any $\theta\in\Theta$, it holds that

$||\Gamma^{n}|| \leq N\max_{1\leq i\leq N}||J_{m_{i}}(\lambda_{i})^{n}||.$

Now, since

$\lim_{n\rightarrow\infty}\sup_{|\lambda_{i}|\leq
1-\alpha}||J_{m_{i}}(\lambda_{i})^{n}||=0$

and $J_{m_{i}}(\lambda_{i})$ is diagonal for $|\lambda_{i}|>1-\alpha$, there
exists a $C\geq 0$ not depending on $\theta$ and $n$ and such that

$||\Gamma^{t}||\leq C|\lambda_{max}(\Gamma)|^{t}$

for all $\theta\in\Theta$. The result then follows from $\sup_{\theta\in
R_{n,0}}|\lambda_{max}(\Gamma)|\leq 1-\log(n)/n$.

For the proof of (e), part (d) and the fact that $\Sigma$ is uniformly bounded
gives us

$\displaystyle\sup_{\theta\in
R_{n,0}}\left|\left|\sum_{t=0}^{n-2}\Gamma^{t}\Sigma\left(\Gamma^{t}\right)^{T}\right|\right|$
$\displaystyle\leq C\sum_{t=0}^{n-2}\sup_{\theta\in
R_{n,0}}||\Gamma^{t}||^{2}\leq C\sum_{t=0}^{n-2}\left(1-\frac{\log
n}{n}\right)^{2t}$ $\displaystyle\leq C\sum_{t=0}^{\infty}\left(1-\frac{\log
n}{n}\right)^{2t}=\frac{C}{(1-(1-\log(n)/n)^{2})}$
$\displaystyle=O\left(\frac{n}{\log n}\right).$

The proof of part (f) is almost the same. Indeed, by the same chain of
inequalities, we find that

$\sup_{\theta\in
R_{n,0}}\left|\left|\sum_{t=0}^{n-2}t\Gamma^{t}\Sigma\left(\Gamma^{t}\right)^{T}\right|\right|\leq\frac{C(1-\log(n)/n)^{2}}{(1-(1-\log(n)/n)^{2})^{2}}=O\left(\frac{n}{\log
n}\right).$

∎

######  Lemma F.2.

We have

$\lim_{n\rightarrow\infty}\sup_{\theta\in
R_{n,0}}\left|\left|H^{-\frac{1}{2}}\left(\sum_{t=0}^{t-2}\Gamma^{t}\Sigma\left(\Gamma^{t}\right)^{T}\right)H^{-\frac{1}{2}}-I\right|\right|=0.$

###### Proof.

Note that it suffices to proof that

$\lim_{n\rightarrow\infty}\sup_{\theta\in
R_{n,0}}\left|\left|\left(\sum_{t=0}^{t-2}\Gamma^{t}\Sigma\left(\Gamma^{t}\right)^{T}\right)^{-\frac{1}{2}}H\left(\sum_{t=0}^{t-2}\Gamma^{t}\Sigma\left(\Gamma^{t}\right)^{T}\right)^{-\frac{1}{2}}-I\right|\right|=0.$

We have

$\displaystyle H$
$\displaystyle=\frac{1}{n}\sum_{t=1}^{n}\sum_{s=1}^{t-1}\Gamma^{t-1-s}\Sigma\left(\Gamma^{t-1-s}\right)^{T}$
$\displaystyle=\frac{1}{n}\sum_{t=0}^{n-2}(n-1-t)\Gamma^{t}\Sigma\left(\Gamma^{t}\right)^{T}$
$\displaystyle=\frac{n-1}{n}\sum_{t=0}^{n-2}\Gamma^{t}\Sigma\left(\Gamma^{t}\right)^{T}-\frac{1}{n}\sum_{t=1}^{n-2}t\Gamma^{t}\Sigma\left(\Gamma^{t}\right)^{T}.$

and, as a result,

$\left(\sum_{t=0}^{t-2}\Gamma^{t}\Sigma\left(\Gamma^{t}\right)^{T}\right)^{-\frac{1}{2}}H\left(\sum_{t=0}^{t-2}\Gamma^{t}\Sigma\left(\Gamma^{t}\right)^{T}\right)^{-\frac{1}{2}}=\frac{n-1}{n}I+M$

where

$M=\left(\sum_{t=0}^{t-2}\Gamma^{t}\Sigma\left(\Gamma^{t}\right)^{T}\right)^{-\frac{1}{2}}\frac{1}{n}\sum_{t=1}^{n-2}t\Gamma^{t}\Sigma\left(\Gamma^{t}\right)^{T}\left(\sum_{t=0}^{t-2}\Gamma^{t}\Sigma\left(\Gamma^{t}\right)^{T}\right)^{-\frac{1}{2}}.$

All that is left to show is therefore that $M$ goes to 0 uniformly over
$R_{n,0}$ for $n$ going to infinity. Similar arguments as in the proof of part
(a) of Lemma F.1 yield

$\liminf_{n\rightarrow\infty}\inf_{\theta\in\Theta}\sigma_{min}\left(\sum_{t=0}^{t-2}\Gamma^{t}\Sigma\left(\Gamma^{t}\right)^{T}\right)>0$

so that

$\sup_{\theta\in\Theta}\left|\left|\left(\sum_{t=0}^{t-2}\Gamma^{t}\Sigma\left(\Gamma^{t}\right)^{T}\right)^{-\frac{1}{2}}\right|\right|<\infty.$

By part (f) of Lemma F.1, we then get

$\sup_{\theta\in R_{n,0}}||M||=o(1).$

∎

Define the diagonal matrix $\Lambda\in\mathbb{R}^{d\times d}$ such that
$\Lambda_{kk}=|\lambda_{i_{k}}|$. The following gives some useful relations
between $I-\Lambda$ and $H$.

######  Lemma F.3.

Let $\Lambda\in\mathbb{R}^{d\times d}$ be as defined above. Then,

(36) $\sup_{\theta\in\Theta}\sup_{t\geq
1}\left|\left|\mathbb{E}\left((I-\Lambda)^{\frac{1}{2}}X_{t,\theta}X_{t,\theta}^{T}(I-\Lambda)^{\frac{1}{2}}\right)\right|\right|<\infty$

and, furthermore,

(37) $\lim_{n\rightarrow\infty}\sup_{\theta\in
R_{n,0}}\left|\left|(I-\Lambda)^{\frac{1}{2}}H(I-\Lambda)^{\frac{1}{2}}-\Sigma_{X}\right|\right|=0$

where
$\textnormal{vec}(\Sigma_{X})=(I-\Lambda)^{\frac{1}{2}}\otimes(I-\Lambda)^{\frac{1}{2}}(I-\Gamma\otimes\Gamma)^{-1}\textnormal{vec}(\Sigma)$
with

$\limsup_{n\rightarrow\infty}\sup_{\theta\in
R_{n,0}}\left\\{\sigma_{min}\left(\Sigma_{X}\right)^{-1}+\sigma_{max}\left(\Sigma_{X}\right)\right\\}<\infty.$

###### Proof.

For any $\theta\in\Theta$ with all diagonal entries of $\Gamma$ having
magnitude greater than $1-\alpha$ and strictly smaller than 1, we have

$\displaystyle\left|\left|\mathbb{E}\left((I-\Lambda)^{\frac{1}{2}}X_{t,\theta}X_{t,\theta}^{T}(I-\Lambda)^{\frac{1}{2}}\right)\right|\right|$
$\displaystyle=\left|\left|(I-\Lambda)^{\frac{1}{2}}\sum_{s=0}^{t-1}\Gamma^{s}\Sigma\left(\Gamma^{s}\right)^{T}(I-\Lambda)^{\frac{1}{2}}\right|\right|$
$\displaystyle\leq||\Sigma||\sum_{s=1}^{t-1}\left|\left|(I-\Lambda)\Gamma^{2s}\right|\right|$
$\displaystyle\leq||\Sigma||\left|\left|(I-\Lambda)(I-\Gamma^{2})^{-1}\right|\right|\left|\left|\left(I-\Gamma^{2t}\right)\right|\right|$
$\displaystyle\leq C\left|\left|(I-\Lambda)(I-\Gamma)^{-1}\right|\right|\leq
Cd$

where

$C=\sup_{\theta\in\Theta}\sup_{t\geq
1}||\Sigma||\left|\left|I+\Gamma\right|\right|\left|\left|\left(I-\Gamma^{2t}\right)\right|\right|<\infty$

and $||(I-\Lambda)(I-\Gamma)^{-1}||\leq d$ follows from Assumption U.5. If
$\Gamma$ is not diagonal the non-diagonal blocks have eigenvalues with
magnitude bounded by $1-\delta$ and hence a similar bound can be achieved. If
$\Gamma$ has a diagonal entry with magnitude $1$, the corresponding diagonal
of $I-\Lambda$ will be 0 and thus the bound above applies for this $\Gamma$
with the corresponding blocks set to 0. This proves (36).

For the proof of (37), simply note that Lemma F.2 and (36) imply that

$\lim_{n\rightarrow\infty}\sup_{\theta\in
R_{n,0}}\left|\left|(I-\Lambda)^{\frac{1}{2}}\left(H-\mathbb{E}(X_{n-1,\theta}X_{n-1,\theta}^{T})\right)(I-\Lambda)^{\frac{1}{2}}\right|\right|=0$

and

$\displaystyle\sup_{\theta\in
R_{n,0}}\left|\left|(I-\Lambda)^{\frac{1}{2}}\mathbb{E}(X_{n-1,\theta}X_{n-1,\theta}^{T})(I-\Lambda)^{\frac{1}{2}}-\Sigma_{X}\right|\right|$
$\displaystyle=$ $\displaystyle\sup_{\theta\in
R_{n,0}}\left|\left|(I-\Lambda)^{\frac{1}{2}}\sum_{s=n-1}^{\infty}\Gamma^{s}\Sigma\left(\Gamma^{s}\right)^{T}(I-\Lambda)^{\frac{1}{2}}\right|\right|\rightarrow
0$

for $n\rightarrow\infty$. It remains to check that $\Sigma_{X}$ is uniformly
bounded and invertible in the limit. Since $\limsup_{n}\sup_{\theta\in
R_{n,0}}\sigma_{max}(\Sigma_{X})<\infty$ follows immediately from (36), we
only need to show the latter. Let $\theta\in R_{n,0}$ be such that $\Gamma$ is
diagonal. Then,

$\displaystyle\sigma_{min}(\Sigma_{X})$
$\displaystyle\geq\sigma_{min}(\Sigma)\sum_{t=1}^{\infty}\sigma_{min}\left((I-\Lambda)\Gamma^{2t}\right)$
$\displaystyle\geq\sigma_{min}(\Sigma)\sum_{t=0}^{\infty}\min_{1\leq k\leq
d}(1-|\lambda_{i_{k}}|)|\lambda_{i_{k}}|^{2t}$
$\displaystyle\geq\sigma_{min}(\Sigma)\sum_{t=0}^{\infty}\frac{\log
n}{n}\left(1-\frac{\log n}{n}\right)^{2t}$
$\displaystyle=\sigma_{min}(\Sigma)\frac{\log n}{n}\left(1-\left(1-\frac{\log
n}{n}\right)^{2}\right)^{-1}$
$\displaystyle\geq\frac{\sigma_{min}(\Sigma)}{2}.$

If $\Gamma$ is non-diagonal, the same bound holds since adding ones on the
super-diagonal does not decrease the minimum singular value. Because $\Sigma$
is uniformly invertible over $\Theta$, the proof is then complete. ∎

The next lemma says that checking whether $A_{n}B_{n}$ converges to the
identity matrix is the same as checking whether $A_{n}B_{n}^{2}A_{n}$
converges to the identity matrix where $A_{n}$ and $B_{n}$ are positive
semidefinite matrices of conforming dimension.

######  Lemma F.4.

Let $\mathcal{I}$ be some index set and consider two families of sequences of
positive semidefinite $d\times d$ matrices,
$(A_{n,i})_{n\in\mathbb{N},i\in\mathcal{I}}$ and
$(B_{n,i})_{n\in\mathbb{N},i\in\mathcal{I}}$. Then, if

$\lim_{n\rightarrow\infty}\sup_{i\in\mathcal{I}}\left|\left|A_{n,i}B_{n,i}^{2}A_{n,i}-I\right|\right|=0,$

it also holds that

$\lim_{n\rightarrow\infty}\sup_{i\in\mathcal{I}}\left|\left|A_{n,i}B_{n,i}-I\right|\right|=0.$

###### Proof.

First, note that

$\liminf_{n\rightarrow\infty}\inf_{i\in\mathcal{I}}\lambda_{min}(A_{n,i}B_{n,i})=c>0.$

Since $A_{n,i}B_{n,i}$ is similar to the positive definite matrix
$A_{n,i}^{\frac{1}{2}}B_{n,i}A_{n,i}^{\frac{1}{2}}$, the contrary would imply
the existence of a sequence $(i_{n})_{n\in\mathbb{N}}\subset\mathcal{I}$ such
that $\sigma_{min}(A_{n,i_{n}}B_{n,i_{n}})\rightarrow 0$ and, thus,
$\sigma_{min}(A_{n,i_{n}}B_{n,i_{n}}^{2}A_{n,i_{n}})\rightarrow 0$ for
$n\rightarrow\infty$ which, of course, is a contradiction.

Now, let $A_{n,i}B_{n,i}=U_{n,i}P_{n,i}$ be a polar decomposition, i.e.,
$U_{n,i}$ is orthogonal and $P_{n,i}$ positive semidefinite. Then,
$P_{n,i}^{2}=A_{n,i}B_{n,i}^{2}A_{n,i}$ which implies that

$\lim_{n\rightarrow\infty}\sup_{i\in\mathcal{I}}\left|\left|P_{n,i}-I\right|\right|=0.$

It therefore suffices to show that $U_{n,i}$ converges to the identity matrix
uniformly over $\mathcal{I}$. Since $U_{n,i}$ is orthogonal, we get

$\left|\left|U_{n,i}-I\right|\right|^{2}=2d-2\text{tr}(U_{n,i})\leq
2d\sup_{1\leq j\leq d}\left|1-\lambda_{j}\left(U_{n,i}\right)\right|.$

Let $U_{n,i}=V_{n,i}D_{n,i}V_{n,i}^{*}$ be an eigendecomposition with
$D_{n,i}$ diagonal and $V_{n,i}$ unitary. Define the Hermitian matrix
$H_{n,i}=V_{n,i}\left(P_{n,i}-I\right)V_{n,i}^{*}$. Denote by $d_{n,i}^{j}$
the $j$’th diagonal of $D_{n,i}$ and by $h_{n,i}^{jk}$ the $jk$’th element of
$H_{n,i}$. Since $H_{n,i}\rightarrow 0$ uniformly over $\mathcal{I}$, for
fixed $\epsilon>0$, we can pick $n_{0}\in\mathbb{N}$ such that

$\sup_{i\in\mathcal{I}}\sup_{1\leq j,k\leq d}|h_{n,i}^{jk}|<\epsilon$

and

$\inf_{i\in\mathcal{I}}\lambda_{min}\left(A_{n,i}B_{n,i}\right)>\frac{c}{2}$

for all $n\geq n_{0}$. We define the complex disk
$D_{r}(x)=\\{y\in\mathbb{C}:|y-x|\leq r\\}$ for any $x\in\mathbb{C}$ and
$r>0$. The matrix $D_{n,i}+D_{n,i}H_{n,i}$ is simlar to $A_{n,i}B_{n,i}$ so,
by the Gershgorin circle theorem, for $1\leq k\leq d$, there exists $1\leq
j\leq d$ such that

$\lambda_{k}\left(A_{n,i}B_{n,i}\right)\in
D_{R_{j}}\left(d_{n,i}^{j}+d_{n,i}^{j}h_{n,i}^{jj}\right)$

where $R_{j}=\sum_{k\neq j}|h_{n,i}^{jk}|$. Recall that $|d_{n,i}^{j}|=1$.
Using the fact that $\lambda_{k}(A_{n,i}B_{n,i})>c/2$ is real and that

$\sup_{i\in\mathcal{I}}\lambda_{k}\left(A_{n,i}B_{n,i}\right)^{2}\leq\sup_{i\in\mathcal{I}}\sigma_{max}\left(A_{n,i}B_{n,i}^{2}A_{n,i}\right)\rightarrow
1$

for $n\rightarrow\infty$, we may therefore assume that $n_{0}$ is large enough
so that

$\sup_{i\in\mathcal{I}}\left|\lambda_{k}\left(A_{n,i}B_{n,i}\right)-1\right|\leq\sup_{i\in\mathcal{I}}\left|d_{n,i}^{j}-\lambda_{k}\left(A_{n,i}B_{n,i}\right)\right|+\epsilon$

for all $n\geq n_{0}$ which implies

$\sup_{i\in\mathcal{I}}\left|d_{n,i}^{j}-1\right|\leq
2\sup_{i\in\mathcal{I}}\left|d_{n,i}^{j}-\lambda_{j}\left(A_{n,i}B_{n,i}\right)\right|+\epsilon\leq
2\sup_{i\in\mathcal{I}}R_{j}+\epsilon+\epsilon\leq 2d\epsilon.$

Thus,

$\sup_{i\in\mathcal{I}}\left|\left|U_{n,i}-I\right|\right|\leq 4d^{2}\epsilon$

for all $n\geq n_{0}$. This completes the proof. ∎

## References

  * [1] Donald WK Andrews  “Exactly median-unbiased estimation of first order autoregressive/unit root models”  In _Econometrica: Journal of the Econometric Society_ JSTOR, 1993, pp. 139–165
  * [2] John Y Campbell and Motohiro Yogo  “Efficient tests of stock return predictability”  In _Journal of financial economics_ 81.1 Elsevier, 2006, pp. 27–60
  * [3] Christophe Cuny, Jérôme Dedecker and Florence Merlevède  “Rates of convergence in invariance principles for random walks on linear groups via martingale methods”  In _Transactions of the American Mathematical Society_ 374.1, 2021, pp. 137–174
  * [4] Robert M De Jong  “Weak laws of large numbers for dependent random variables”  In _Annales d’Économie et de Statistique_ JSTOR, 1998, pp. 209–225
  * [5] Juan J Dolado and Helmut Lütkepohl  “Making Wald tests work for cointegrated VAR systems”  In _Econometric reviews_ 15.4 Taylor & Francis, 1996, pp. 369–386
  * [6] James Duffy and Jerome Simons  “Cointegration without Unit Roots”, 2023
  * [7] Graham Elliott  “On the robustness of cointegration methods when regressors almost have unit roots”  In _Econometrica_ 66.1 JSTOR, 1998, pp. 149–158
  * [8] Graham Elliott and James H Stock  “Inference in time series regression when the order of integration of a regressor is unknown”  In _Econometric theory_ 10.3-4 Cambridge University Press, 1994, pp. 672–700
  * [9] Massimo Franchi and Søren Johansen  “Improved inference on cointegrating vectors in the presence of a near unit root using adjusted quantiles”  In _Econometrics_ 5.2 MDPI, 2017, pp. 25
  * [10] Bruce E Hansen  “The grid bootstrap and the autoregressive model”  In _Review of Economics and Statistics_ 81.4 MIT Press 238 Main St., Suite 500, Cambridge, MA 02142-1046, USA journals …, 1999, pp. 594–607
  * [11] Michael Jansson and Marcelo J Moreira  “Optimal inference in regression models with nearly integrated regressors”  In _Econometrica_ 74.3 Wiley Online Library, 2006, pp. 681–714
  * [12] Hiroaki Kaido, Francesca Molinari and Jörg Stoye  “Confidence intervals for projections of partially identified parameters”  In _Econometrica_ 87.4 Wiley Online Library, 2019, pp. 1397–1432
  * [13] Hiroaki Kaido, Francesca Molinari, Jörg Stoye and Matthew Thirkettle  “Calibrated Projection in MATLAB: Users’ Manual”  In _arXiv preprint arXiv:1710.09707_ , 2017
  * [14] Maximilian Kasy  “Uniformity and the delta method”  In _Journal of Econometric Methods_ 8.1 De Gruyter, 2019
  * [15] Alexandros Kostakis, Tassos Magdalinos and Michalis P Stamatogiannis  “Robust econometric inference for stock return predictability”  In _The Review of Financial Studies_ 28.5 Oxford University Press, 2015, pp. 1506–1553
  * [16] Anton Rask Lundborg, Rajen D Shah and Jonas Peters  “Conditional independence testing in hilbert spaces with applications to functional data analysis”  In _arXiv preprint arXiv:2101.07108_ , 2021
  * [17] Tassos Magdalinos and Peter Charles Bonest Phillips  “Econometric inference in matrix vicinities of unity and stationarity”  In _Working Paper_ , 2020
  * [18] Adam McCloskey  “Bonferroni-based size-correction for nonstandard testing problems”  In _Journal of Econometrics_ 200.1 Elsevier, 2017, pp. 17–35
  * [19] Anna Mikusheva  “Uniform inference in autoregressive models”  In _Econometrica_ 75.5 Wiley Online Library, 2007, pp. 1411–1452
  * [20] Peter CB Phillips  “On confidence intervals for autoregressive roots and predictive regression”  In _Econometrica_ 82.3 Wiley Online Library, 2014, pp. 1177–1195
  * [21] Peter CB Phillips and Tassos Magdalinos  “Econometric inference in the vicinity of unity”  In _Singapore Management University, CoFie Working Paper_ 7, 2009
  * [22] Peter CB Phillips and Tassos Magdalinos  “Limit theory for moderate deviations from a unit root”  In _Journal of Econometrics_ 136.1 Elsevier, 2007, pp. 115–130
  * [23] D. Pollard  “Convergence of Stochastic Processes”, Springer Series in Statistics  Springer New York, 1984
  * [24] Hiro Y Toda and Taku Yamamoto  “Statistical inference in vector autoregressions with possibly integrated processes”  In _Journal of econometrics_ 66.1-2 Elsevier, 1995, pp. 225–250