# LidarAugment: Searching for Scalable 3D LiDAR Data Augmentations

Zhaoqi Leng1∗, Guowang Li1, Chenxi Liu1, Ekin Dogus Cubuk2, Pei Sun1, Tong
He1,
Dragomir Anguelov1 and Mingxing Tan1 1 Waymo Research, 2 Google Brain, ∗
<EMAIL_ADDRESS>
###### Abstract

Data augmentations are important in training high-performance 3D object
detectors for point clouds. Despite recent efforts on designing new data
augmentations, perhaps surprisingly, most state-of-the-art 3D detectors only
use a few simple data augmentations. In particular, different from 2D image
data augmentations, 3D data augmentations need to account for different
representations of input data and require being customized for different
models, which introduces significant overhead. In this paper, we resort to a
search-based approach, and propose _LidarAugment_ , a practical and effective
data augmentation strategy for 3D object detection. Unlike previous approaches
where all augmentation policies are tuned in an exponentially large search
space, we propose to factorize and align the search space of each data
augmentation, which cuts down the 20+ hyperparameters to 2, and significantly
reduces the search complexity. We show LidarAugment can be customized for
different model architectures with different input representations by a simple
2D grid search, and consistently improve both convolution-based
UPillars/StarNet/RSN and transformer-based SWFormer. Furthermore, LidarAugment
mitigates overfitting and allows us to scale up 3D detectors to much larger
capacity. In particular, by combining with latest 3D detectors, our
LidarAugment achieves a new state-of-the-art 74.8 mAPH L2 on Waymo Open
Dataset.

##  I Introduction

Data augmentations are widely used in training deep neural networks. In
particular, for autonomous driving, many data augmentations are developed to
improve data efficiency and model generalization. However, most recent 3D
object detectors only use a few basic data augmentation operations such as
rotation, flip and ground-truth sampling [1, 2, 3, 4, 5, 6, 7]. This is in a
surprising contrast to 2D image recognition and detection, where much more
sophisticated 2D data augmentations are commonly used in modern image-based
models [8, 9, 10, 11, 12, 13]. In this paper, we aim to answer: _is it
practical to adopt more advanced 3D data augmentations to improve modern 3D
object detectors, especially for high-capacity models?_

Figure 1: Model scaling with LidarAugment on Waymo Open Dataset. Baseline
augmentations are from the prior art of [14]. When scaling up UPillars to
UPillars-L, our LidarAugment improves both models, and the gains are more
significant for the larger model, thanks to its customizable regularization.
More results in Table IV. Figure 2: Visualizing LidarAugment. (a) all data
augmentation operations used in LidarAugment. For non-global operations, we
highlight the augmented parts in red (boxes). (b) occlusion introduced by data
augmentation, e.g., paste a car object, is handled by removing overlapping
rays in range view based on distance. We show point clouds and the
corresponding range images with (bottom)/without (top) removing overlapping
rays in the range view.

The main challenge of adopting advanced 3D data augmentations is that 3D
augmentations are often sensitive to input representations and model capacity.
For example, range image based models and point cloud based models require
different types of data augmentation due to different input representations.
High capacity 3D detectors are typically prone to overfitting and require
stronger overall data augmentation compared to lite models with fewer
parameters. Therefore, tailoring each 3D augmentation for different models is
necessary. However, the search space scales exponentially with respect to the
number of hyperparameters, which leads to significant search cost. Recent
studies [15, 16] attempt to address these challenges by using efficient search
algorithms. Those approaches typically construct a fixed search space, and run
a complex search algorithms (such as population-based search [17]) to find a
data augmentation strategy for a model. However, our studies reveal that the
search spaces used in prior works are suboptimal. Despite having complex
search algorithms, without a systematic way to define a good search space, we
cannot unleash the potential of a model.

In this paper, we propose LidarAugment, a simplified search-based approach for
3D data augmentations. Unlike previous methods that rely on complex search
algorithms to explore an exponentially large search space, our approach aims
to define a simplified search space that contains a variety of data
augmentations but has minimal (i.e. two) hyperparameters, such that users can
easily customize a diverse set of 3D data augmentations for different models.

Specifically, we construct the LidarAugment search space by first factorizing
a large search space based on operations and exploring each sub search space
with a per-operation search. Then, we normalize and align the sub search space
for each data augmentation to form the LidarAugment search space. The final
LidarAugment search space contains only two shared hyperparameters:
$m\in[0,\infty)$ controls the normalized magnitude and $p\in[0,1]$ controls
the probability of applying each data augmentation policies. Our LidarAugment
search space significantly simplifies prior works [15] by cutting down the
number of hyperparameters to two, a 15$\times$ reduction in number of
hyperparameters.

Despite only having two hyperparamters, our LidarAugment search space contains
a variety of existing 3D data augmentations, such as drop/paste 3D bounding
boxes, rotate/scale/dropping points, and copy-paste objects and backgrounds.
In addition, LidarAugment supports coherent augmentation across both point and
range view representations, which generalizes to multi-view 3D detectors.

We perform extensive experiments on the Waymo Open Dataset [18] and
demonstrate LidarAugment is effective and generalizes well to different model
architectures (convolutions-based and transformer-based), different input
views (3D point view and range image), and different temporal scales (single
and multi frames). Notably, LidarAugment advances state-of-the-art (SOTA)
transformer-based SWFormer by 1.4 mAPH on the test set. Furthermore,
LidarAugment provides customizable regularization, which allows us to scale up
3D object detectors to much higher capacity without overfitting. As summarized
in Figure 1, LidarAugment consistently improves UPillars models, and the
performance gains are particularly large for high-capacity models. Our
contributions can be summarized as:

  1. 1.

New insight: we reveal that common 3D data augmentation search spaces are
suboptimal and should be tailored for different models.

  2. 2.

LidarAugment: we propose the LidarAugment search space, which supports jointly
optimizing 10 augmentation policies with only two hyperparameters (15$\times$
reduction compares to prior works), offering diverse yet practical
augmentations. In addition, we develop a new method to coherently augment both
point and range-view input representations.

  3. 3.

State-of-the-art performance: LidarAugment consistently improves both
convolution-based UPillars/StarNet/RSN and attention-based SWFormer. With
LidarAugment, we achieve new state-of-the-art results on Waymo Open Dataset.
In addition, LidarAugment enables model scaling to achieve much better quality
for high-capacity 3D detectors.

##  II Related works

Data augmentation. Data augmentation is widely used in training deep neural
networks In particular, for 3D object detection from point clouds, several
global and local data augmentations, such as rotation, flip, pasting objects,
and frustum noise, are used to improve model performance [19, 1, 20, 2, 4, 21,
15, 22, 23, 24]. However, as 3D data augmentations are sensitive to model
architectures and capacity, it often requires extensive manual tuning to use
these augmentations. Therefore, most existing 3D object detectors [2, 6, 25,
26, 14] only adopt a few simple augmentations, such as flip and shift pixels.

Several recent works attempt to use range images for multi-view 3D detection,
but very few augmentations are developed for range images. [5] attempts to
paste objects in the range image without handling occlusions. Our Paste Box
augmentation support coherently augmenting both range-view and point-view
input data while handling occluded objects in a simple way (more details in
Figure 2), which enables more realistic augmented scenes and enriches the data
augmentations for multi-view 3D detectors.

Learning data augmentation policies. Designing good data augmentation normally
requires manual tuning and domain expertise. Several search-based approaches
have been proposed for 2D images, such as AutoAugment [9], RandAugment [12],
and Fast AutoAugment [27]. Our LidarAugment is inspired by RandAugment in the
sense that we both try to construct a simplified search space. However, unlike
2D image augmentations, where a search space works well for many models, we
reveal that existing search space for 3D detection tasks are suboptimal, which
motivates us to propose the first systematical method to define search spaces
for 3D detection tasks.

On the other hand, for 3D detection, PPBA [15] and PointAugment [16] propose
efficient learning-based data augmentation frameworks for 3D point clouds.
However, both works require users to run a complex algorithm on an
exponentially large but not well-designed search space. In contrast, our work
provides a systematical framework to design a simple and more effective search
spaces with only two hyperparameters.

##  III LidarAugment

In this section, we first introduce data augmentation policies used in
LidarAugment. Next, we analyze the performance of each data augmentation
policy on Waymo Open Dataset [18]. Finally, we propose a systematic approach
to progressively design 3D augmentation search space.

###  III-A Data augmentations for point clouds and range images.

3D point cloud and 2D range image are two different representations of LiDAR
data. Despite being the native representation of LiDAR data, data
augmentations for range image is not well studied compared to point clouds.
Here, we revisit data augmentations for point clouds, and introduce a new
method for coherently applying data augmentation to both point clouds and
range images.

Augmenting point clouds. We follow the implementation of data augmentation
policies described in recent studies [1, 15, 28], which contain global
operations (rotate, scale, translate, flip, and drop points) and local
operations (drop boxes, paste boxes, swap background, drop points and add
feature noise in a frustum), shown in Figure 2 (a).

Augmenting range images. Different from sparse 3D point representation, pixels
in range image are compact. Data augmentations, such as pasting objects and
swap background, disturb the compact structure of range representation. Here,
we propose a novel approach to coherently augment both 3D point view and 2D
range view by leveraging the bijective property between point clouds and range
images, while account for occlusion.

First, we transform the range image pixels to point cloud based on $(x,y,z)$
coordinates. To preserve the bijective mapping between a pixel in a range
image and a point in the corresponding point clouds, we concatenate the (row,
column) index of each pixel in the range image as additional features before
scattering pixels to 3D. After performing data augmentation in the point
representation, we transform the augmented point clouds back to the range view
by scattering each point to a pixel in a 2D image based on its (row, column)
index.

Leveraging the compactness of range images. Coherently augmenting both range
and point views leads to more realistic augmented scenes. Because each pixel
in a range image corresponds to a unique ray from LiDAR, overlapping pixels in
the range view represent that the same light ray penetrates trough multiple
surfaces. When this happens, we compare the distance among overlapping pixels
in the range view and keep the pixel that is closest to the ego vehicle. This
effectively removes occluded points in both the range and point views, as
shown in Figure 2 (b).

Policy |  Hyperparameters | WOD (Veh./Ped.) | mAP L1
---|---|---|---
No Aug |  - | - | 60.2
Drop Box |  Probability | $p/p$ | 66.0 (+5.8)
Number of boxes | $2m/2.8m$
Paste Box |  Probability | $1.4p/p$ | 66.6 (+6.4)
Number of boxes | $3.2m/4.4m$
Swap Background |  Probability | $0.6p$ | 63.6 (+3.4)
Global Rot |  Probability | $1.4p$ | 73.3 (+13.1)
Max rotation angle | $0.22\pi m$
Global Scale |  Probability | $p$ | 66.0 (+5.8)
Scaling factor | $0.036m$
Global Drop |  Probability | $p$ | 64.9 (+4.7)
Drop ratio | $1-0.18m$
Frustum Drop |  Probability | $p$ | 64.1 (+3.9)
Theta angle width | $0.1\pi m$
Phi angle width | $0.1\pi m$
R distance | $75-7.5m$
Drop ratio | $1-0.1m$
Frustum Noise |  Probability | $0.6p$ | 65.1 (+4.9)
Theta angle width | $0.14\pi m$
Phi angle width | $0.14\pi m$
R distance | $75-10.5m$
Max noise level | $0.14m$
Global Translate |  Probability | $1.4p$ | 67.5 (+7.3)
Stdev. of noise (x, y) | $0.66m$
Global Flip |  Probability | $p$ | 69.0 (+8.8)

TABLE I: Aligned search spaces and performance. The search space of each
hyperparameter for Waymo Open Dataset (WOD) for UPillars is listed. $(p,m)$
are two global hyperparameters to control all data augmentation policies.
After aliging the search space, the optimal $(p,m)$ for each data augmentation
are $(0.5,5)$. The probability of each policy is clipped to [0, 1]. The min R
distance is clipped to 0. The maximum rotation angle is clipped to [0, $\pi$].
The maximum flip probability is clipped to 0.5. The ratio of dropped points
are clipped to [0, 0.8]. The theta angle and phi angle are clipped to [0,
$\pi$] and [0, $2\pi$], respectively.  Figure 3: UPillars architecture. Input
points are processed by two full connected layers with channel size (MLP0,
MLP1) before voxelized into pillars. The bird-eye-view pillars are processed
by a Res U-Net, where the channel size and number of blocks at each resolution
are (Ci, Bi). CenterNet detection head, box regression, and other attributes
regression heads are applied to the output of U-Net. For both models, we use
the same voxel size $0.32m$ and range $81.92m$.

###  III-B Effects of each data augmentation.

In this section, we assess the effects of each data augmentation policy on
Waymo Open Dataset [18]. To benchmark the policies, we develop a UPillars
architecture, which is based on the popular PointPillars [2], but incorporates
recent optimizations in architecture design, i.e., unet backbone [29], and
center net detection head [30].

Datasets and training. Waymo Open Dataset [18] contains 798 and 202 training
and validation sequences. For the following studies, we train UPillars with
batch size 64, Adam optimizer [31] and cosine decay learning rate with max
learning rate 3e-3 and total step 80000.

Effect of each data augmentations. We factorize the LidarAugment search space
into per-policy sub search space and show the UPillars performance when
trained using only one policy on Waymo Open Dataset (WOD) in Table I.
Interestingly, on WOD, the most effective data augmentation technique is
global rotation, whereas on KITTI [32], pasting ground truth bounding boxes is
commonly regarded as the most effective data augmentation [1]. A closer look
at the statistics of the two datasets reveals that, on average, each KITTI
LiDAR frame contains about five objects, whereas, each frame in WOD on average
contains more than 50 objects. Thus, pasting ground truth objects has larger
impact on KITTI, due to the significantly lower object density, than WOD. On
the other hand, smaller global rotation angle $\pi/4$ is commonly used when
training KITTI dataset, but we find much stronger rotation $\pi$ is preferred
for WOD.

###  III-C Defining LidarAugment search space.

As indicated from the previous section, different from RandAugment for 2D
images [12], naively using the same search space across different datasets is
suboptimal, which is a unique challenge for 3D detection tasks. To mitigate
this new challenge, we propose to factorize the whole search space and align
each data augmentation based on its optimal hyperparameters.

Align the search space. Using global hyperparameters to control all data
augmentations requires normalizing the search domain of each hyperparameter.
Without normalization, the same global magnitude could lead to an aggressive
application of one data augmentation, and an insufficient application of
another data augmentation. To align the search domain of each data
augmentation policy, we train a UPillars model on a given dataset while only
applying a single data augmentation policy at each time. Based on the optimal
values of hyperparameters, we rescale the search domain of each hyperparameter
in a data augmentation policy from $[0,arbitrary\\_value]$ to
$[0,optimal\\_value]$ such that the optimal value for each hyperparameter
corresponds to the same global magnitude or probability hyperparameter. Since
each data augmentation policy contains multiple parameters, to save cost, we
perform a small scale 2D grid search to scale the probability and magnitudes
of all hyperparameters in each sub search space.

⬇

augmentations = [

DropBox, PasteBox, SwapBackground, GlobalRot,

GlobalScale, GlobalDrop, FrustumDrop ,

FrustumNoise, GlobalTranslate, GlobalFlip]

def lidaraugment(m, p, input_frame):

for aug in augmentations:

aug.set_magnitude(m)

aug.set_probability(p)

input_frame = aug.transform(input_frame)

return input_frame

Figure 4: Pseudo Python code for LidarAugment.

Here, we use Global Translate as an example. We define the initial search
domain for the probability of applying Global Translate data augmentation to
be $\\{0.3,0.5,0.7,0.9\\}$, and the domain for the magnitude of translation
noise to be $\\{0.9,1.5,2.1,2.7,3.3,3.9\\}$. If the optimal values
$(p_{\text{noise}},m_{\text{noise}})=(0.7,3.3)$, we rescale the search domain
to $(p_{\text{noise}},m_{\text{noise}})=(1.4p,0.66m)$ such that when the
global hyperparameters $(p,m)=(0.5,5)$, the hyperparameters for Global
Translate are optimal. Details about all the hyperparameters are listed in
Table I. The LidarAugment pseudocode is shown in Figure 4.

##  IV Experiments

In this section, we first introduce our experimental setups. Then we show
LidarAugment significantly improves the performance for both convolution-based
and attention-based models. Lastly, we show the model scaling results,
followed by ablations studied on different models and datasets.

###  IV-A Experimental setup

Our experiments are mostly based on Waymo Open Dataset [18] (WOD) where the
main metric is mAPH L2, with additional ablation studies on nuScene [33]. We
evaluate LidarAugment on a variety of 3D object detectors, as well as
different model sizes. For fair comparison, we strictly follow the original
training settings for each model, and only replace the baseline augmentation
with our new LidarAugment. We train UPillars using Adam optimizer [31] and
apply cosine learning rate with max learning rate 1e-3, total 16e4 steps and
batch size 64.

Method | Type | mAPH | Vehicle AP/APH 3D | Pedestrian AP/APH 3D
---|---|---|---|---
L2 | L1 | L2 | L1 | L2
P.Pillars [2] † | conv | 51.9 | 63.3/62.7 | 55.2/54.7 | 68.9/56.6 | 60.4/49.1
CenterPoint [25] | conv | 67.1 | 76.6/76.1 | 68.9/68.4 | 79.0/73.4 | 71.0/65.8
RSN_3f [6] | conv | 68.1 | 78.4/78.1 | 69.5/69.1 | 79.4/76.2 | 69.9/67.0
PVRCNN++ [34] | conv | 69.1 | 79.3/78.8 | 70.6/70.2 | 81.8/76.3 | 73.2/68.0
UPillars-L† | conv | 60.0 | 69.5/69.0 | 61.5/61.0 | 70.4/66.1 | 63.0/59.0
UPillars-L(+LA) | conv | 71.0 | 79.5/79.0 | 71.9/71.5 | 81.5/77.3 | 74.5/70.5
SST_1f [26] | attn | 63.4 | 74.2/73.8 | 65.5/65.1 | 78.7/69.6 | 70.0/61.7
SST_3f [26] | attn | 69.5 | 77.0/76.6 | 68.5/68.1 | 82.4/78.0 | 75.1/70.9
SWFormer [14] | attn | 70.9 | 79.4/78.9 | 71.1/70.6 | 82.9/79.0 | 74.8/71.1
SWFormer(+LA) | attn | 72.8 | 80.9/80.4 | 72.8/72.4 | 84.4/80.7 | 76.8/73.2

TABLE II: WOD validation-set results. _LA_ denotes our LidarAugment, _conv_
denotes convolutional networks, and _attn_ denotes attention-based transformer
models. LidarAugment improves both types of models, and achieves the best
results among each category. † model is trained using augmentations shown in
prior art [14].

###  IV-B LidarAugment achieves new state-of-the-art results

Table II compares the validation set results on Waymo Open dataset. Our
LidarAugment significantly improves both convolution-based and transformer-
based models. In particular, by scaling up the basic UPillars, our
LidarAugment achieves 71.0 mAPH of L2 on UPillars-L, which is 1.9 AP better
than the previous best convolution-based 3D detector PVRCNN++ [34]. Notably,
the latest transformer-based SWFormer [14] already uses 4 strong data
augmentation policies, i.e., rotation (probability 0.74, yam angle uniformly
sampled from [$-\pi$, $\pi$]), random flip (probability 0.5), randomly scaling
the world (scaling factor uniformly sampled from $[0.95,1.05]$, and randomly
drop points (drop probability 0.05), where the rotation angle and flip
probability are maxed out. Despite that, LidarAugment still outperforms
SWFormer by 1.9 AP, establishing a new state-of-the-art result for single-
modal models without ensemble or test time augmentation on Waymo Open Dataset.

Table III compares the test-set results among latest models. Compared to the
latest SWFormer, our LidarAugment improves the test-set L2 mAPH by 1.4 AP,
outperforming all prior arts by a large margin.

Method | mAPH | Vehicle AP/APH 3D | Pedestrian AP/APH 3D
---|---|---|---
L2 | L1 | L2 | L1 | L2
P.Pillars [2] † | 55.1 | 68.6/68.1 | 60.5/60.1 | 68.0/55.5 | 61.4/50.1
CenterPoint [25] | 69.1 | 80.2/79.7 | 72.2/71.8 | 78.3/72.1 | 72.2/66.4
RSN_3f [6] | 69.7 | 80.7/80.3 | 71.9/71.6 | 78.9/75.6 | 70.7/67.8
PVRCNN++ [34] | 71.2 | 81.6/81.2 | 73.9/73.5 | 80.4/75.0 | 74.1/69.0
SST_TS_3f [26] | 72.9 | 81.0/80.6 | 73.1/72.7 | 83.1/79.4 | 76.7/73.1
SWFormer [14] | 73.4 | 82.9/82.5 | 75.0/74.7 | 82.1/78.1 | 75.9/72.1
SWFormer(+LA) | 74.8 | 84.0/83.6 | 76.3/76.0 | 83.1/79.3 | 77.2/73.5

TABLE III: WOD test-set results. LidarAugment (LA) significantly improves
detection performance for SWFormer, and achieves new state-of-the-art mAPH L2.

###  IV-C LidarAugment enables better model scaling

Scaling up model capacity is a common approach to achieve better performance,
but large 3D object detectors often suffer from overfitting. Table IV shows
scaling results, where UPillars-L is a larger model with more layers and
channels than UPillars, detailed in Figure 3.

Here, we adopt the strong data augmentations used in the latest SWFormer as
Baseline (see subsection IV-B). As shown in Table IV, with Baseline
augmentations, UPillars-L does not benefit much from its significantly larger
capacity. In fact, several metrics, such as Veh/Ped L1 AP, even become worse
(e.g. 69.3/70.3 for UPillar-L vs 72.1/72.3 for UPillars). We observe the
training loss of UPillars-L is much smaller compared to loss of UPillars,
indicating severe overfitting.

On the other hand, LidarAugment achieves much better performance on larger
models, especially on the most challenging metric, i.e., +7.3AP for 3D L2 mAPH
as shown in Figure 1. Perhaps surprisingly, although the baseline UPillar
(mAPH=57.8) is much worse than latest 3D detectors, the final performance of
UPillar-L (+ LidarAugment) is actually competitive with the latest SWFormers,
i.e., their mAPH are 71.0 vs. 72.8. This opens up new research opportunities
on exploring much larger and higher performance 3D detectors in the future.

| Veh/Ped AP L1 | Veh/Ped APH L2
---|---|---
UPillars | UPillars-L | UPillars | UPillars-L
BaseAugment | 72.1/72.3 | 69.3/70.3 | 63.5/52.1 | 61.0/59.0
LidarAugment | 77.1/77.5 | 79.5/81.6 | 68.5/58.9 | 71.5/70.5

TABLE IV: UPillars scaling results on WOD.

###  IV-D LidarAugment supports different representations

Different from 2D image models, 3D detectors are more diverse and could
utilize different input representations due to the additional dimensionality
and sparsity of point cloud data. Other than UPillars and SWFormer, which are
both pillar-based architectures and taking 3D sparse points as inputs, we
further demonstrate LidarAugment generalizes to other input representations.
First, StarNet [35] is a point-based detector which directly processes raw
points in 3D to detect objects. RSN, on the other hand, utilize multi-view
property of point clouds and takes both range images and 3D sparse points as
inputs. However, due to the lack of multi-view data augmentations in prior
works, RSN only utilize two simple augmentations, i.e. random flip and
rotation.

Model | Augmentation | Vehicle | Pedestrian
---|---|---|---
StarNet [35] | baseline | 58.2 | 71.9
+LidarAugment | 61.6 | 74.2
RSN-1frame [6] | baseline | 75.2 | 77.2
+LidarAugment | 75.8 | 79.0
RSN-3frame [6] | baseline | 77.0 | 79.1
+LidarAugment | 77.7 | 80.6
SWFormer [14] | baseline | 79.4 | 82.9
+LidarAugment | 80.9 | 84.4

TABLE V: LidarAugment improves various models. Startnet is a point-based
detector. RSN is a range image and pillar-based detector. Results are WOD L1
AP.

LidarAugment is a general method which supports augmenting different views of
point clouds, including range images, as explained in subsection III-A. Table
V shows the performance of LidarAugment on point-based StarNet, range image
based RSN, and transformer-based SWFormer. In general, our LidarAugment
improves all kinds of 3D detectors, sometimes by a large margin.

###  IV-E Abletion studies: comparing to other approaches.

In this section, we show LidarAugment outperforms other common data
augmentation approaches on UPillars.

Manually tuned data augmentation. Due to the complexity of search space scales
exponentially with respect to the number of parameters, commonly used data
augmentation strategies often consists of few data augmentation operations.
Here, we benchmark two sets of data augmentation strategies used in training
high-performance 3D detectors. First, we adopt random flip (probability 0.5)
and rotation (probability 0.5, yaw angle uniformly sampled from [$-\pi/4$,
$\pi/4$]) data augmentations used in training RSN [6]. Then we benchmark more
advanced and stronger data augmentation strategy used in training SWFormer
[14], detailed in subsection IV-B. Our results show both data augmentation
strategies significantly improved UPillars performances, about +10 AP for
Vehicle and Pedestrian 3D L1 AP, when compared to the no augmentation
baseline, shown in Table VI. However, tuning the data augmentation
hyperparameters is challenging, e.g., if we only search 4 values for each
hyperparameter, the number of searches of 5 hyperparameters exceeds 1000.

UPillars (AP Level 1) | Hparams | Vehicle | Pedestrian
---|---|---|---
No Augmentation | - | 58.0 | 62.4
Rotate & Flip [6] | 3 | 70.8 | 70.0
Rotate, Flip, Scale, Drop points [14] | 5 | 72.1 | 72.3
PPBA [15] | 29 | 71.6 | 72.6
LidarAugment | 2 |  77.1 (+5.0) |  77.5 (+4.9)

TABLE VI: LidarAugment outperforms common data augmentation strategies.
UPillars L1 APs on Waymo Open Dataset validation set are reported.
LidarAugment requires the least number of hyperparameters (Hparams) but
achieves the best results compared to manually designed and automl-based data
augmentations strategies.

AutoML-based data augmentation. To alleviate the challenge of exponentially
large search space, population-based training is proposed to tune
hyperparameters in data augmentations online [17, 11, 15]. We follow the
implementation of progressive-population based data augmentation (PPBA) [15]
and use the same sets of data augmentation policies and search space. We set
population size 16, generation step 4000, perturbation and exploration rate to
0.2. Our results, in Table VI, show PPBA significant outperforms the no
augmentation baseline. Despite PPBA introduces significantly more data
augmentation policies, it is on par with manually tuned data augmentations,
which only contains 4 policies.

We find the search space of PPBA is suboptimal after inspecting the search
domain of each hyperparameter. For example, the maximum rotation angle for
global rotation in the PPBA search space is $\pi/4$, a common value used for
KITTI dataset. However, $\pi/4$ is insufficient compared to the tailored max
rotation angle $\pi$ used in our LidarAugment. Surprisingly, a single well-
tuned global rotation augmentation achieves L1 mAP 73.3, shown in Table I,
which outperforms PPBA with L1 mAP 72.1 over vehicle and pedestrian tasks.
Although PPBA algorithm is more efficient than grid search and contains
diverse augmentation policies, the suboptimal search domain of rotation angle
restricts the performance of PPBA, which highlights the importance of
tailoring 3D detection search space.

LidarAugment Alternatively, LidarAugment mitigates both the curse of
dimensionality and suboptimal search domain issues by aligning and scaling the
magnitude and probability of each data augmentation policy. This significantly
reduces the search complexity (only 2 hyperparameters) while allowing
exploration of a larger hyperparameter space. As indicate in Table VI,
LidarAugment significantly outperforms both manually designed and AutoML-based
data augmentation strategies by about 5 AP for both vehicle and pedestrian
detection tasks and only requires a simple grid search of two hyperparameters.

###  IV-F Generalize to nuScenes dataset

To further validate our method, we evaluate LidarAugment on a different
dataset: nuScenes [33]. For simplicity, we adopt the same training settings as
Waymo Open Dataset, but reduce the voxel size to 0.25 and the total training
steps by half for faster training. We use the same baseline augmentation as
SWFormer, and redefine LidarAugment search space for nuScenes following
subsection III-C. Table VII shows LidarAugment is a general approach, which
outperforms the baseline augmentation by a large margin on nuScenes.

UPillars | mAP | NDS
---|---|---
Rotate, Flip, Scale, Drop points [14] | 40.6 | 48.2
LidarAugment |  46.7 (+6.1) |  53.4 (+5.2)

TABLE VII: nuScenes validation-set results.

##  V Conclusion

In this paper, we propose _LidarAugment_ , a scalable and effective 3D
augmentation approach for 3D object detection. Based on the insight that 3D
data augmentations are sensitive to model architecture and capacity, we
propose a simplified search space, which contains two hyperparameters to
control a diverse set of augmentations. LidarAugment outperforms both manually
tuned and existing search-based data augmentation strategies by a large
margin. Extensive studies show that LidarAugment generalizes to convolution
and attention-based architectures, as well as point-based and range-based
input representations. More importantly, LidarAugment significantly simplifies
the search process for 3D data augmentations and opens up exciting new
research opportunities, such as model scaling in 3D detection. With
LidarAugment, we demonstrate new state-of-the-art 3D detection results on the
challenging Waymo Open Dataset.

## References

  * [1] Y. Yan, Y. Mao, and B. Li, “Second: Sparsely embedded convolutional detection,” _Sensors_ , vol. 18, no. 10, p. 3337, 2018.
  * [2] A. H. Lang, S. Vora, H. Caesar, L. Zhou, J. Yang, and O. Beijbom, “Pointpillars: Fast encoders for object detection from point clouds,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 12 697–12 705.
  * [3] Y. Zhou, P. Sun, Y. Zhang, D. Anguelov, J. Gao, T. Ouyang, J. Guo, J. Ngiam, and V. Vasudevan, “End-to-end multi-view fusion for 3d object detection in lidar point clouds,” in _Conference on Robot Learning_. PMLR, 2020, pp. 923–932.
  * [4] S. Shi, C. Guo, L. Jiang, Z. Wang, J. Shi, X. Wang, and H. Li, “Pv-rcnn: Point-voxel feature set abstraction for 3d object detection,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2020, pp. 10 529–10 538.
  * [5] Z. Liang, M. Zhang, Z. Zhang, X. Zhao, and S. Pu, “Rangercnn: Towards fast and accurate 3d object detection with range image representation,” _arXiv preprint arXiv:2009.00206_ , 2020.
  * [6] P. Sun, W. Wang, Y. Chai, G. Elsayed, A. Bewley, X. Zhang, C. Sminchisescu, and D. Anguelov, “Rsn: Range sparse net for efficient, accurate lidar 3d object detection,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2021, pp. 5725–5734.
  * [7] A. Bewley, P. Sun, T. Mensink, D. Anguelov, and C. Sminchisescu, “Range conditioned dilated convolutions for scale invariant 3d object detection,” in _Conference on Robot Learning_. PMLR, 2021, pp. 627–641.
  * [8] H. Zhang, M. Cisse, Y. N. Dauphin, and D. Lopez-Paz, “mixup: Beyond empirical risk minimization,” _arXiv preprint arXiv:1710.09412_ , 2017.
  * [9] E. D. Cubuk, B. Zoph, D. Mane, V. Vasudevan, and Q. V. Le, “Autoaugment: Learning augmentation policies from data,” _arXiv preprint arXiv:1805.09501_ , 2018.
  * [10] S. Yun, D. Han, S. J. Oh, S. Chun, J. Choe, and Y. Yoo, “Cutmix: Regularization strategy to train strong classifiers with localizable features,” in _Proceedings of the IEEE/CVF international conference on computer vision_ , 2019, pp. 6023–6032.
  * [11] D. Ho, E. Liang, X. Chen, I. Stoica, and P. Abbeel, “Population based augmentation: Efficient learning of augmentation policy schedules,” in _International Conference on Machine Learning_. PMLR, 2019, pp. 2731–2741.
  * [12] E. D. Cubuk, B. Zoph, J. Shlens, and Q. V. Le, “Randaugment: Practical automated data augmentation with a reduced search space,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops_ , 2020, pp. 702–703.
  * [13] Z. Zhong, L. Zheng, G. Kang, S. Li, and Y. Yang, “Random erasing data augmentation,” in _Proceedings of the AAAI conference on artificial intelligence_ , vol. 34, no. 07, 2020, pp. 13 001–13 008.
  * [14] P. Sun, M. Tan, W. Wang, C. Liu, F. Xia, Z. Leng, and D. Anguelov, “Swformer: Sparse window transformer for 3d object detection in point clouds,” _European Conference on Computer Vision (ECCV)_ , 2022.
  * [15] S. Cheng, Z. Leng, E. D. Cubuk, B. Zoph, C. Bai, J. Ngiam, Y. Song, B. Caine, V. Vasudevan, C. Li _et al._ , “Improving 3d object detection through progressive population based augmentation,” in _European Conference on Computer Vision_. Springer, 2020, pp. 279–294.
  * [16] R. Li, X. Li, P.-A. Heng, and C.-W. Fu, “Pointaugment: an auto-augmentation framework for point cloud classification,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2020, pp. 6378–6387.
  * [17] M. Jaderberg, V. Dalibard, S. Osindero, W. M. Czarnecki, J. Donahue, A. Razavi, O. Vinyals, T. Green, I. Dunning, K. Simonyan _et al._ , “Population based training of neural networks,” _arXiv preprint arXiv:1711.09846_ , 2017\.
  * [18] P. Sun, H. Kretzschmar, X. Dotiwalla, A. Chouard, V. Patnaik, P. Tsui, J. Guo, Y. Zhou, Y. Chai, B. Caine _et al._ , “Scalability in perception for autonomous driving: Waymo open dataset,” in _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , 2020, pp. 2446–2454.
  * [19] X. Chen, H. Ma, J. Wan, B. Li, and T. Xia, “Multi-view 3d object detection network for autonomous driving,” in _Proceedings of the IEEE conference on Computer Vision and Pattern Recognition_ , 2017, pp. 1907–1915.
  * [20] B. Yang, W. Luo, and R. Urtasun, “Pixor: Real-time 3d object detection from point clouds,” in _Proceedings of the IEEE conference on Computer Vision and Pattern Recognition_ , 2018, pp. 7652–7660.
  * [21] Y. Chen, V. T. Hu, E. Gavves, T. Mensink, P. Mettes, P. Yang, and C. G. Snoek, “Pointmixup: Augmentation for point clouds,” in _European Conference on Computer Vision_. Springer, 2020, pp. 330–345.
  * [22] J. S. Hu and S. L. Waslander, “Pattern-aware data augmentation for lidar 3d object detection,” in _2021 IEEE International Intelligent Transportation Systems Conference (ITSC)_. IEEE, 2021, pp. 2703–2710.
  * [23] J. Choi, Y. Song, and N. Kwak, “Part-aware data augmentation for 3d object detection in point cloud,” in _2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2021, pp. 3391–3397.
  * [24] M. Reuse, M. Simon, and B. Sick, “About the ambiguity of data augmentation for 3d object detection in autonomous driving,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2021, pp. 979–987.
  * [25] T. Yin, X. Zhou, and P. Krahenbuhl, “Center-based 3d object detection and tracking,” in _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , 2021, pp. 11 784–11 793.
  * [26] L. Fan, Z. Pang, T. Zhang, Y.-X. Wang, H. Zhao, F. Wang, N. Wang, and Z. Zhang, “Embracing single stride 3d object detector with sparse transformer,” _arXiv preprint arXiv:2112.06375_ , 2021.
  * [27] S. Lim, I. Kim, T. Kim, C. Kim, and S. Kim, “Fast autoaugment,” _Advances in Neural Information Processing Systems_ , vol. 32, 2019.
  * [28] Z. Leng, S. Cheng, B. Caine, W. Wang, X. Zhang, S. Jonathon, M. Tan, and A. Dragomir, “Pseudoaugment: Learning to use unlabeled data for data augmentation in point clouds,” in _European Conference on Computer Vision_. Springer, 2022, pp. 279–294.
  * [29] O. Ronneberger, P. Fischer, and T. Brox, “U-net: Convolutional networks for biomedical image segmentation,” in _International Conference on Medical image computing and computer-assisted intervention_. Springer, 2015, pp. 234–241.
  * [30] T. Yin, X. Zhou, and P. Krahenbuhl, “Center-based 3d object detection and tracking,” in _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , 2021, pp. 11 784–11 793.
  * [31] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.
  * [32] A. Geiger, P. Lenz, C. Stiller, and R. Urtasun, “Vision meets robotics: The kitti dataset,” _The International Journal of Robotics Research_ , vol. 32, no. 11, pp. 1231–1237, 2013.
  * [33] H. Caesar, V. Bankiti, A. H. Lang, S. Vora, V. E. Liong, Q. Xu, A. Krishnan, Y. Pan, G. Baldan, and O. Beijbom, “nuscenes: A multimodal dataset for autonomous driving,” in _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , 2020, pp. 11 621–11 631.
  * [34] S. Shi, L. Jiang, J. Deng, Z. Wang, C. Guo, J. Shi, X. Wang, and H. Li, “Pv-rcnn++: Point-voxel feature set abstraction with local vector representation for 3d object detection,” _arXiv preprint arXiv:2102.00463_ , 2021.
  * [35] J. Ngiam, B. Caine, W. Han, B. Yang, Y. Chai, P. Sun, Y. Zhou, X. Yi, O. Alsharif, P. Nguyen _et al._ , “Starnet: Targeted computation for object detection in point clouds,” _arXiv preprint arXiv:1908.11069_ , 2019\.


# Identification testing via sample splitting – an application to Structural
VAR models

Katarzyna Maciejowska<EMAIL_ADDRESS>Department of
Operations Research and Business Intelligence, Faculty of Management, Wrocław
University of Science and Technology, 50-370 Wrocław, Poland

(This version: )

###### Abstract

In this article, a novel identification test is proposed, which can be applied
to parameteric models such as Mixture of Normal (MN) distributions, Markow
Switching(MS), or Structural Autoregressive (SVAR) models. In the approach, it
is assumed that model parameters are identified under the null whereas under
the alternative they are not identified. Thanks to the setting, the Maximum
Likelihood (ML) estimator preserves its properties under the null hypothesis.
The proposed test is based on a comparison of two consistent estimators based
on independent subsamples of the data set. A Wald type statistic is proposed
which has a typical $\chi^{2}$ distribution. Finally, the method is adjusted
to test if the heteroscedasticity assumption is sufficient to identify
parameters of SVAR model. Its properties are evaluated with a Monte Carlo
experiment, which allows non Gaussian distribution of errors and mis-specified
VAR order. They indicate that the test has an asymptotically correct size.
Moreover, outcomes show that the power of the test makes it suitable for
empirical applications.

###### keywords:

identification test, Maximum Likelihood, Wald test, Structural VAR,
heterosecedasticity

††journal: Renewable Energy

##  1 Introduction

Identifications is one of the crucial assumptions, which needs to be fulfilled
in order to use properly the Maximum Likelihood estimation method (see Hsiao
(1983), Rothenberg (1971), among others). For majority of econometric models,
this property is easily proved and hence is not a concern for researchers.
Unfortunately, there are families of models, in which the identifications is
not guarantee, for example Structural Vector Autoregressive (SVAR) models
(Sims (1980), Lütkepohl (2005), Kikian and Lütkepohl (2017)), nonlinear regime
switching models such as Smooth Transition Regression (STR) or mixture models
(McLachlan et al. (2018), van Dijk et al. (2002)).

In Rothenberg (1971), two types of identification types are discussed: global
and local. Both definitions are based on a parametric likelihood function of a
random variable, $y\in R^{M}$. Lets denote by $l(\theta,y)$ a log-likelihood
functions computed for the parameter vector $\theta$. Then it is said that the
the true parameter vector $\theta_{0}$ is globally identifiable if for

$\forall_{\theta\in\Theta}\forall_{y\in
R^{K}}:\theta_{0}=\theta\Leftrightarrow l(\theta_{0};y)=l(\theta;y),$ (1)

Global identifications implies that the ML estimator is unique and hence under
regularity conditions, the Maximum Likelihood (ML) estimator of $\theta$ is
consistent. Although desirable, the global identification is not necessary to
estimate the model parameters. The sufficient condition is local
identification. The parameter vector $\theta_{0}$ is locally identifiable if

$\forall_{\varepsilon}\forall_{\|\theta_{0}-\theta\|<\varepsilon}\forall_{y\in
R^{K}}:\theta_{0}=\theta\Leftrightarrow l(\theta_{0},y)=l(\theta,y).$ (2)

When the local identification is ensured, the parameter space, $\Theta$, can
be restricted in such way that the parameters becomes unique. One of the
examples of models which are identifiable only localy is a Mixture of Normal
distributions (MN) model. It can be shown that when the ordering of components
in the mixture is changed, which is called a ’label switching’, the parameter
vectors change as well but the value of the likelihood functions remain
unaffected. This problem can be solved by assuming that the probability of the
first state is larger than 0.5 or that the mean of the first component is
larger than the second one.

Unfortunately, there are cases, in which the local identification is not
satisfied. For example, when the MN model is fitted to data stemming from a
simple normal distribution or when the SVAR model is applied without imposing
sufficient restrictions. In such case, there are infinitely many parameter
vectors, which provide the same value of the log-likelihood function. As the
result, it is not guaranteed that the ML estimator converge to $\theta_{0}$ as
the sample size increases. Moreover, as shown by Rothenberg (1971), the
information matrix for $\theta_{0}$ is singular and hence the typical ML test
(Wald, LM or LR) either cannot be computed or have an unknown limiting
distribution.

In the literature, there are a few identification tests dedicated to
particular econometric problems. For example, in case of switching regression
such as Markov Switching (MS), STR or structural break modes, it is known as a
linearity test because for linear processes the models are not identified.
Luukkonen et al. (1988) use a Taylor expansion of the STR model to verify if
the underlying process is linear. Andrews (1993) explores _Sup_ transformation
of classical ML statistics and proves its asymptotic properties. Hansen (1992)
proposed a test for MS models, which is based on LR statistic but has a
nonstandard distribution. For SVAR models, the problem is directly called
identification and has been studied by Lütkepohl and Woźniak (2021) and
Lütkepohl et al. (2021). However, there is no single testing approach, which
will be suitable to more general family of parametric models and would be
straightforward to apply. Andrews and Ploberger (1994) and Hansen (1996)
develop tests for cases when the nuisance parameters are presented only under
the alternative hypothesis. Unfortunately, their statistics have nonstandard
asymptotic distributions, which cannot be tabulated accept special cases.
Hence, applying these tests requires conducting additional, complex
simulations. Moreover and cannot be applied to some models, for example SVAR.

In this research, a novel approach is proposed for testing identification,
which can be applied to different parametric models. It differs from the
existing methods, for example a recent test of Lütkepohl et al. (2021) or
previous Hansen (1996) and Andrews and Ploberger (1994) tests, in terms of the
hypothesis setting. In the presented approach, similar to Rossi and Inoue
(2011), it is assumed that the model is identified under the null and not-
identified under the alternative. Hence the test supports identification
hypothesis, as long as there are no strong evidence against it. This is
convenient for researchers who want to use for MN or SVAR models when the data
does not reject them. Additionally, it can be notice that under the null,
there should be not estimation problems and the information matrix resulting
from ML method is non-singular.

The testing procedure explores the fact that under the null any consistent
estimator should converge to its true value, which is uniquely identified. In
order to obtain two independent estimators, the sample is split randomly into
two disjoint sets. Next the two estimators are compared with the Wald type of
test. This approach has a few advantages. First, it can be applied to any
parametric model if sufficient data is available. Second, since the Wald type
of test is based on the information matrix rather than its inverse, the
statistic could be computed even when the Hessian is singular (under the
alternative). Finally, the test has a well known $\chi^{2}$ distribution and
hence is does not require any additional simulations and is straightforward to
apply.

This article demonstrates how the approach can be adopted for testing
identifications in SVAR models. Recent publications (see Lanne and Lütkepohl
(2008), Lanne et al. (2010) and Lütkepohl and Natšunajev (2017), among others)
show that statistical data properties, such as heteroscedasticity, can be
explored to identify structural parameters without imposing additional
restriction. The crucial assumption in this procedure is that shock variances
change in different proportions. Although this assumption is very important,
it cannot be directly tested. Some papers still use the standard Wald or
Likelihood Ratio (LR) tests to evaluate the model, for example Lanne et al.
(2010), Lütkepohl and Velinov (2021). Unfortunately, a $\chi^{2}$ distribution
does not provide a proper approximation of their asymptotic behavior. In other
research (Lütkepohl and Woźniak (2021)), Bayesian methods are used to verify
the model assumptions. Recently, Lütkepohl et al. (2021) proposed a direct
frequentist test for assessing the change of volatility in SVAR models. This
article explores the properties of eigen values of covariance matrices to
construct a test with an asymptotic $\chi^{2}$ distribution. Albeit its
novelty and straightforward application, the new test is designed only for
models with discrete shift of variances and it is not examined how it works
for Markov-switching (MS) or STR types of heteroscedasticity. Moreover, as
under the null the model is not identified, it is not clear how to state the
hypothesis when the evaluated process has a dimension larger than two. The
authors propose to apply sequential testing. A method proposed in this article
does not suffer from the above problems. There is one null hypothesis, under
which the structure of SVAR models is well identified. The test statistic has
an asymptotic $\chi^{2}$ distribution and can be directly applied to different
types of heteroscedasticity.

Finally, the performance of a newly proposed identification test is evaluated
with a MC experiment. In order to provide results comparable with the previous
outcomes of Lütkepohl et al. (2021), the same data generating processes are
considered. They cover the cases of non-Gaussian distributions, a misspecified
order of VAR model and processes of a dimension larger than two. The results
indicate that the proposed test has a proper size. Its power increases as the
p-values form different sample splitting are combined via a harmonic mean.
Moreover, it is robust to model misspecification and works well for higher
dimension processes. Finally, the approach is applied to verify the
identification of Kilian (2009). The results confirm findings of Lütkepohl and
Natšunajev (2014) and indicate that there are no strong evidences that the
model is not well identified via heteroscedasticity.

The article is structured as follows. In the first Section 2 the proposed
identification test is described. Next, its application to SVAR model is
described in Section 3. The performance of the test is evaluated by a Monte
Carlo experiment, which results are shown in Section 4. Section 5 presents an
exmaple of an empirical application of the testing approach and Section 6
concludes.

##  2 Identification test

Lets consider a parametric model, which is described by a log-likelihood
function $L(\theta,Y)=\sum_{t=1}^{T}l(\theta,y_{t})$, where
$Y=\\{y_{1},...,y_{T}\\}$ is the information set and $T$ is the sample size.
The parameter vector $\theta\in R^{M}$. In the test, two hypothesis are
considered: under the null the model parameters are assumed to be identifiable
whereas under the alternative they lack identification. As the result, under
$H_{0}$ (if the likelihood functions satisfied standard regularity conditions
(Greene (2018))) the ML estimator becomes consistent and asymptotically
normal. These properties are not maintained for $H_{1}$, under which the
estimator does not converge in probability to $\theta_{0}$. The proposed
approach explores the consistency property to derive the test statistic.

In order to build the test statistic, the sample is split into two disjoint
parts and the estimators coming from these two sub-samples are compared with
the Wald type of test. The testing procedure consists of the following steps

  1. 1.

Randomly split the information set into two disjoint parts of equal sizes:
$Y_{1}$, $Y_{2}$, such that $Y_{1}\cup Y_{2}=Y$, $Y_{1}\cap Y_{2}=\varnothing$

  2. 2.

Compute ML estimator for the whole information set and the two parts:
$\hat{\theta}$, $\hat{\theta}_{1}$, $\hat{\theta}_{2}$, respectively

  3. 3.

Compute the average Hessian of the full sample log-likelihood functions for
$\hat{\theta}$

$H_{T}(\hat{\theta})=\frac{1}{T}\sum_{t=1}^{T}{\frac{\partial^{2}l(\theta;y_{t})}{\partial\theta\partial\theta^{\prime}}(\hat{\theta})}$
(3)

  4. 4.

Compute the Wald type statistic

$W=T/2(\hat{\theta}_{1}-\hat{\theta}_{2})^{\prime}(-H_{T}(\hat{\theta})/2)(\hat{\theta}_{1}-\hat{\theta}_{2})$
(4)

It can be shown that under the null the statistic $W$ converges to the
$\chi^{2}(M)$ distribution, where $M$ is the size of the parameter vector
$\theta$.

######  Theorem 1.

Suppose a model described by a log-likelihood function $l(\theta,y_{t})$
satisfies regularity conditions, which ensure consistency of the ML estimator
and its asymptotic normality. Then under the null hypothesis of parameter
identification, the test statistic

$W=T/2(\hat{\theta}_{1}-\hat{\theta}_{2})^{\prime}(-H_{T}(\hat{\theta})/2)(\hat{\theta}_{1}-\hat{\theta}_{2})$

has an asymptotic $\chi^{2}(M)$ distribution.

###### Proof.

First, it can be observed that $\hat{\theta}_{1}$ and $\hat{\theta}_{2}$ are
independent because they come from two different, randomly drown sub-samples.
Second, as the model is identifiable under the null then the estimators are
consistent and the CLT implies that

$\sqrt{T/2}(\hat{\theta}_{i}-\theta_{0})\rightarrow_{d}N(0,(-H(\theta_{0}))^{-1}),$
(5)

where $i=1,2$ and $H(\theta_{0})=p\lim_{T\to\infty}H_{T}(\theta_{0})$. The
difference of two independent normal random variables is also normal, with the
variance being the sum of the individual variances. Hence

$\sqrt{T/2}(\hat{\theta}_{1}-\theta_{1})\rightarrow_{d}N(0,2(-H(\theta_{0}))^{-1}),$
(6)

Finally, as $\hat{\theta}$ is a consistent ML estimator with
$\hat{\theta}\rightarrow_{p}\theta_{0}$ then also
$H_{T}(\hat{\theta})\rightarrow_{p}H_{T}(\theta_{0})$ and finally

$H_{T}(\hat{\theta})\rightarrow_{p}H(\theta_{0}).$ (7)

As a result,

$T/2(\hat{\theta}_{1}-\hat{\theta}_{2})^{\prime}(-H_{T}(\hat{\theta})/2)(\hat{\theta}_{1}-\hat{\theta}_{2})\rightarrow_{d}\chi^{2}(M).$
(8)

∎

The statistic $W$ has an intuitive interpretation. It shows that if the model
is identifiable, the estimators stemming from two datasets should be similar
and hence their difference converges to zero. It can be noticed that if the
parameter estimates are not identifiable, the two estimators
$\hat{\theta}_{1}$ and $\hat{\theta}_{2}$ do not necessary converge to the
same value and hence their difference preserves not negligible even for large
samples.

####  2.0.1 Behavior of test under the alternative

Behavior of the test under the alternative is a complex and nontrivial
problem. When ML estimator is not consistent then the limit of $\hat{\theta}$
does not exist, and so the limits of $\hat{\theta}_{1}$ and
$\hat{\theta}_{2}$. Then the difference $\hat{\theta}_{1}-\hat{\theta}_{2}$
does not converge in probability to a zero vector.

######  Theorem 2.

Suppose there exists a set $A\subset R^{M}$, which includes at least two
elements, such that if and only if $\tilde{\theta}\in A$, the log-likelihood
$l(\tilde{\theta};y)=l(\theta_{0};y)$. So the parameters in $A$ are observable
equivalent to the true parameter vector $\theta_{0}$. Then the difference of
two independent ML estimators, $\hat{\theta}_{1}-\hat{\theta}_{2}$, does not
converge in probability to zero.

###### Proof.

First, lets denote by $e$, a distance of a parameters vector $\theta$ to a set
$A$ defined as

$e=\theta-arg\min_{\tilde{\theta}\in A}{|\theta-\tilde{\theta}|}$

Then the parameter can be expressed as a sum $\theta=\tilde{\theta}+e$, with
$\tilde{\theta}\in A$. Since $\theta_{0}\in A$ then the ML estimator will
converge to $A$, meaning that $p\lim\hat{e}=0$.

Lets now consider the difference $\hat{\theta}_{1}-\hat{\theta}_{2}$. It could
be expressed as

$\hat{\theta}_{1}-\hat{\theta}_{2}=\tilde{\theta}_{1}-\tilde{\theta}_{2}+\hat{e}_{1}-\hat{e}_{2}.$

Since $p\lim\hat{e}=0$ then also $p\lim(\hat{e}_{1}-\hat{e}_{2})=0$. So the
difference $\hat{\theta}_{1}-\hat{\theta}_{2}$ will converge in probability to
zero if and only if $p\lim(\tilde{\theta}_{1}-\tilde{\theta}_{2})=0$. The
parameters $\tilde{\theta}_{1}$ and $\tilde{\theta}_{2}$ are equally likely,
hence

$prob(|\tilde{\theta}_{1}-\tilde{\theta}_{2}|>\varepsilon)=\frac{||\tilde{\theta}_{1}-\tilde{\theta}_{2}|>\varepsilon|}{|A|}.$

It can be noticed that since $A$ contains at least two elements and
$\tilde{\theta}_{1}$ and $\tilde{\theta}_{2}$ are independent then for
$\varepsilon=0$ the probability
$prob(\tilde{\theta}_{1}=\tilde{\theta}_{2})=0$ and hence

$prob(|\tilde{\theta}_{1}-\tilde{\theta}_{2}|>0)\geq 0.5.$

Hence there exists $\varepsilon$ such that

$prob(|\tilde{\theta}_{1}-\tilde{\theta}_{2}|>\varepsilon)\nrightarrow 0.$

This implies that the difference $\hat{\theta}_{1}-\hat{\theta}_{2}$ does not
converge in probability to zero. ∎

Although the distance between the two estimators $\hat{\theta}_{1}$ and
$\hat{\theta}_{2}$ does not converge to zero, it does not directly imply that
the $W$ statistic under the alternative grows big enough to exceed the
critical values. This is due to the fact that the lack of identification is
associated with a singularity of a second derivative matrix, $H$ Rothenberg
(1971). Therefore, the behavior of the test statistic should be examined
separately for different types of models and identification problems.

###  2.1 Multiple splitting

In the proposed approach, the data split is random – not driven by the data
itself. Hence it seems natural to repeat the splitting for w few times. Lets
denote by $N$ the number of times the data is divided and the statistic
$W_{n}$ with a corresponding p-value, $v_{n}$ is calculated. Here
$n=1,2,...,N$.

Following Vovk and Wang (2020), the p-values will be combined via averaging.
As stated in the paper, the harmonic mean is a robust approach, which yields
good results for both correlated and uncorrelated inputs. The main concern
associated with application of the theory is the fact that the $p$-values are
obtained with the asymptotic rather then actual distribution of the test
statistic. Hence, some of the assumptions stated by Vovk and Wang (2020) are
not met. In particular, for small samples the $p$-values are not uniformly
distributed. Therefore, the harmonic mean is applied here to the truncated
sample $v_{n}\in V$, which includes only the observations of $v_{n}$ staying
between the 20th and 80th percentiles of its empirical distribution. The
truncation removes from the average elements, which are either very big (when
$v_{n}$ is close to zero), or close to one. This change may have a positive
impact on both the size and the power of the test, as it reduces symmetrically
the input noise.

As the result, the averaged $p$-value, $\bar{v}_{N}$, is calculated according
to the following formula

$\bar{v}_{N}=a_{N}\left(\frac{1}{\tilde{N}}\sum_{v_{n}\in
V}{\frac{1}{v_{n}}}\right)^{-1}.$ (9)

It can be noticed that once the set $V$ is changed, the inverse of the
$p$-values is averaged over a different number of observations, $\tilde{N}$.
Here, $\tilde{N}=0.6N$. Finally, as stated by Vovk and Wang (2020), the
harmonic mean should be scaled by a factor $a_{N}$, which is

$a_{N}=\frac{(z^{*}+N)^{2}}{(z^{*}+1)N,}$ (10)

where $z^{*}$ is the unique solution ($z\in(0,\infty)$) to the equation

$z^{2}=N((z+1)ln(z+1)-z).$

The scaling factor is shown to $a_{N}\leq e\ln N$ and $a_{N}/\ln{N}\to 1$.
Unfortunately, the convergence is very slow hence $a_{N}$ is approximated
numerically separately for each $N$. The table with values of $a_{N}$ for
selected number of averaged $p$-values can be found in Vovk and Wang (2020).

It is expected that the theoretical properties of the truncated harmonic mean
can be different from the untruncated one. However simulation results indicate
its validity and importance for small samples, in which the distribution of
$p$-values deviates significantly from the uniform one. The robustness of the
averaging method to the selection of $V$ is discussed and evaluated in Section
4.4.

##  3 Test for an identification of a SVAR model

###  3.1 SVAR model: identification and estimation

The identification problem is one of the central issues in SVAR modeling. The
Vector Autoregressive (VAR) model describes the join behavior of a set of time
series. Lets define a vector of endogenous variables as $Y_{t}\in R^{K}$. The
reduced form of SVAR model takes the following form

$Y_{t}=\mu+\sum_{p=1}^{P}A_{p}Y_{t-p}+\varepsilon_{t},$ (11)

where $\mu$ is a $K\times 1$ vector of intercepts, $A_{p}$ are $K\times K$
matrices describing the autoregresive structure and $\varepsilon_{t}$ is a
$K\times 1$ vector of residuals. Since $\varepsilon_{t}$ are forecast errors,
they are allowed to be correlated and hence their variance, $\Sigma$, is not
diagonal.

The reduced form model, although useful for modeling of the expected value of
endogenous variables, is not suited for risk analysis – for example impulse
responses or simulations. Therefor, a structural extension of the model (11)
is used. The SVAR model may be described in various ways (see Lütkepohl (2005)
for an extension discussion). Here, the B-model is presented, in which the
residuals, $\varepsilon_{t}$, are a linear function of structural shocks,
$u_{t}$. The structural shocks are assumed to be uncorrelated and have a
diagonal variance-covariance matrix, $\Lambda$. Then

$\varepsilon_{t}=Bu_{t}.$ (12)

The matrix $h$ is called _an instantaneous effect matrix_ , because it
describes how structural shocks affect endogenous variables in the current
time period, $t$. For example, $B_{ij}$ describes the impact of the $j$-th
structural shock, $u_{j,t}$ on the $i$-th element of $Y_{t}$. Notice that
equation (12) implies that

$\Sigma=B\Lambda B^{\prime},$ (13)

so there is a direct relationship between $B$ and the variance of errors
$\varepsilon_{t}$. It is typically assumed that either structural shocks have
an identity variance-covariance matrix, $\Lambda=I$, or the diagonal elements
of $B$ are equal to one. Here, the first approach is adopted.

Unfortunately, the parameters of SVAR model cannot be directly estimated
because there are infinitely many matrices $B$, which satisfies (13). As shown
in the literature (see Lütkepohl, 2005, for a comprehensive discussion on VAR
models), the structural model requires estimation of $K^{2}$ elements of the
$B$ matrix. At the same time, the variance-covariance matrix $\Sigma$ of the
reduced form consists of only $K(K+1)/2$ distinct parameters. Therefore
additional $K(K-1)/2$ restrictions need to be imposed, which will restrict the
parameter space.

The identification restrictions may stem either from the theoretical knowledge
of the analyzed economic problem or from statistical properties of the data.
Lanne and Lütkepohl (2008) shows that if there is a shift in the variance of
structural shocks, one can recover the $B$ matrix without imposing additional
constraints. Lanne and Lütkepohl (2008) assumes that the shift is due to a
structural change, which appears at time $T_{c}$, and hence for $t\leq T_{c}$
the variance of shocks $\Sigma_{u}=I$, whereas for $t>T_{c}$ there is
$\Sigma_{u}=\Lambda$. They stated that if the diagonal elements of $\Lambda$
are all different then the model parameters are identifiable. Unfortunately,
this assumption could not be tested, because if some of the elements were
identical, the model would become not identifiable and standard ML based
statistics could not be used.

Using the results from previous Section 2, a two step procedure, similar to
Lütkepohl et al. (2021), is proposed to estimate SVAR model. In the first
step, the parameters of the reduced form of the model (11) are estimated with
the Generalized Least Square (GLS) method. As the result, the residuals
$\hat{\varepsilon}_{t}=Y_{t}-\hat{\mu}-\sum_{p=1}^{P}\hat{A}_{p}Y_{t-p}$ are
computed. In the second step, the structure model of errors,
$\hat{\varepsilon}_{t}$, is estimated with the ML method and its
identification is tested. In ML, it is assumed that shocks are independent and
Gaussian. As the result, the residuals of the reduce form VAR follows:

$\varepsilon_{t}\sim\left\\{\begin{array}[]{lcl}N(0,BB^{\prime})&\mbox{for}&t\leq
T_{c}\\\ N(0,B\Lambda B^{\prime})&\mbox{for}&t>T_{c}.\\\ \end{array}\right.$
(14)

The log-likelihood function used in ML becomes

$l(\theta)=-T\log\det(B)-\frac{T-T_{c}}{2}\log\det(\Lambda)-\frac{1}{2}\sum_{t=1}^{T_{c}}\varepsilon_{t}^{\prime}(BB^{\prime})^{-1}\varepsilon_{t}-\frac{1}{2}\sum_{t=T_{c}+1}^{T}\varepsilon_{t}^{\prime}(B\Lambda
B^{\prime})^{-1}\varepsilon_{t},$ (15)

where $\theta=(vec(B)^{\prime},\lambda_{1},...,\lambda_{K})^{\prime}$ is a
vector of parameters (’vec’ is an operator, which stacks columns of the given
matrix into a single vector). As shown by Lütkepohl et al. (2021), if we
denote $\Sigma_{1}=BB^{\prime}$ and $\Sigma_{2}=B\Lambda B^{\prime}$ then
diagonal elements of $\Lambda$ are the eigenvalues of the matrix
$\Sigma_{2}\Sigma_{1}^{-1}$. The columns of the matrix $B=[b_{1},...,b_{K}]$
are proportional to the eigenvectors satisfying the equation:

$\left(\Sigma_{1}^{-1}-\lambda_{k}\Sigma^{-1}_{2}\right)b_{k}=0,$ (16)

for $k=1,2,..,K$. As shown by Lanne et al. (2010), when the eigenvalues are
all distinct then the matrix $B$ is unique up to a permutation and a sign of
its columns. Therefore, in order to ensure that the two estimators:
$\hat{\theta}_{1}$ and $\hat{\theta}_{2}$ corresponds to each other, it is
assumed that $\lambda$’s are in a descending order
($\lambda_{1}\geq\lambda_{2}\geq...\geq\lambda_{K}$) and that an element of
each column $b_{k}$ with the largest absolute value is positive.

###  3.2 Identification test

In order apply the testing procedure described in Section, the following pair
of hypothesis is considered

$H_{0}:\forall k\in\\{1,...,K-1\\}:\lambda_{k}\neq\lambda_{k+1}$
$H_{1}:\exists k\in\\{1,...,K-1\\}:\lambda_{k}=\lambda_{k+1}$

It can be noticed that the hypothesis are opposite to those stated by
Lütkepohl et al. (2021), which assumes that the model is identifiable under
the alternative. As the result, Lütkepohl et al. (2021) needed to consider a
set of different null hypothesis for models with $K>2$ and test the
identifications sequentially. Here, the null hypothesis is unique and hence it
is sufficient to compute the test only once.

The testing procedure consists of the following steps

  1. 1.

Split randomly the sample $S=\\{1,2,...,T\\}$ into $S_{1}$, $S_{2}$ in such
way that: $|S_{1}|=|S_{2}|=T/2$ and $|S_{1}\leq T_{c}|=|S_{2}\leq
T_{c}|=T_{c}/2$, where $|S_{i}\leq T_{c}|$ is the number of observations such
that $S_{i}\leq T_{c}$.

  2. 2.

Estimate the model parameters for the whole sample, $\hat{\theta}$, and the
two sub-samples $\hat{\theta}_{1}$, $\hat{\theta}_{2}$.

  3. 3.

Compute the average second derivative of the Log Likelihood function for
$\hat{\theta}$: $H_{T}(\hat{\theta})$

  4. 4.

Compute the test statistic

$W=-T/4(\hat{\theta}_{1}-\hat{\theta}_{2})^{\prime}H_{T}(\hat{\theta})(\hat{\theta}_{1}-\hat{\theta}_{2})$
(17)

and compare it with $\chi^{2}(K^{2}+K)$. When the null is not rejected, there
is no evidence that the model is not identifiable.

When the underlying distribution of structural shocks is non Gaussian then the
likelihood function is misspecified. In particular, $H_{T}(\hat{\theta})$ may
not converge to its true values. For heavy tails distributions, such as
$t$-Student distribution, the parameter variances can be underestimated and
the $W$ statistic may reject the null too often. Since the Hessian of
$l(\theta)$ depends on the inverse of the square of error variance, whereas
the covariance of parameters depend on the fourth central moment of errors,
the discrepancy between $E\varepsilon_{t}^{4}$ and
$(E\varepsilon_{t}^{2})^{2}$ is partially responsible for the misperformance
of the test. Therefor, a kurtosis – $\kappa$ – of the underlying DGP seems
crucial for adjusting the test to the non-Gaussian distributions. Here, a
modified version of the $W$ statistic is proposed, in which

$\tilde{H}_{T}(\theta,\kappa)=H_{T}(\theta)/(\kappa/3).$ (18)

and

$W(\hat{\kappa})=-T/4(\hat{\theta}_{1}-\hat{\theta}_{2})^{\prime}\tilde{H}_{T}(\hat{\theta},\hat{\kappa})(\hat{\theta}_{1}-\hat{\theta}_{2})=-T/4(\hat{\theta}_{1}-\hat{\theta}_{2})^{\prime}H_{T}(\hat{\theta})(\hat{\theta}_{1}-\hat{\theta}_{2})/(\hat{\kappa}/3),$
(19)

where $\hat{\kappa}$ is a consistent estimator of $\kappa$. In this article,
an estimator discussed by Schott (2001) and applied in Lütkepohl et al. (2021)
is used in order to make the results comparable. It could be noticed that for
normally distributed residuals, $\kappa=3$ and $W=W(3)$. Moreover, for any
consistent estimator of kurtosis, such that $\hat{\kappa}\to_{p}\kappa$,
$W(\hat{\kappa})\to_{d}\chi^{2}(K^{2}+K)$.

Finally, the steps(1)-(4) can be repeated for $N$ times, as discussed in
Section 2, and the average $p$-value, $\bar{v}_{N}$ can be computed according
to (9). In may happen that $\bar{v}_{N}>1$. In such case, the $\bar{v}_{N}$ is
set to equal to one.

##  4 Small sample properties of the identification test

In order to assess the finite sample properties of the proposed test, a Monte
Carlo (MC) experiment is conducted. It is assumed that structural shocks,
$u_{t}$ are i.i.d. white noise processes, with a covariance matrix
characterized by a structural change in period $T_{c}=\tau T$. Similar to
Lütkepohl et al. (2021), bivariate and three-dimensional GDPs are considered:

  * 1.

$\mathbf{DGP1}$: In the model $y_{t}=\varepsilon_{t}$ and
$\varepsilon_{t}=Bu_{t}$. The structural shocks $u_{t}\in R^{2}$ are
independent Gaussian or $t(5)$ distributions with the covariance matrix
$\Lambda=diag(\lambda_{1},\lambda_{2})$ and the instantaneous effect matrix is
$B=I_{2}$. The structural change occurs at $\tau=0.5$. Analogous to Lütkepohl
et al. (2021), two values of $\Lambda$ are considered:
$(\lambda_{1},\lambda_{2})=(2,2),(2,1)$.

  * 2.

$\mathbf{DGP2}$: A two-dimensional VAR(2) model

$y_{t}=\left[\begin{array}[]{c}0.190\\\ 0.523\\\
\end{array}\right]+\left[\begin{array}[]{cc}0.-0.036&-0.705\\\ -0.093&1.211\\\
\end{array}\right]y_{t-1}+\left[\begin{array}[]{cc}0.090&0.796\\\
-0.085&-0.276\\\ \end{array}\right]y_{t-2}+\varepsilon_{t}.$

The residuals $\varepsilon_{t}=Bu_{t}$. The structural shocks $u_{t}\in R^{2}$
are independent Gaussian with $\Lambda=diag(\lambda_{1},\lambda_{2})$ and
$(\lambda_{1},\lambda_{2})=(0.5,0.5),(0.5,0.1)$. The structural change occurs
at $\tau=0.3$ and the instantaneous effect matrix is

$B=\left[\begin{array}[]{cc}0.317&1.059\\\ 0.242&-0.450\\\
\end{array}\right].$

  * 3.

$\mathbf{GDP3}$: A three dimensional VAR(0) model with $y_{t}=\varepsilon_{t}$
and $\varepsilon_{t}=Bu_{t}$. The structural shocks $u_{t}\in R^{3}$ are
independent Gaussian with $\Lambda=diag(\lambda_{1},\lambda_{2},\lambda_{3})$
and $(\lambda_{1},\lambda_{2},\lambda_{3})=(3,2,1),(3,2,2),(2,2,2)$. The
structural change occurs at $\tau=0.5$ and the instantaneous effect matrix
$B=I_{3}$.

  * 4.

$\mathbf{GDP4}$: A three dimensional VAR(0) model with $y_{t}=\varepsilon_{t}$
and $\varepsilon_{t}=Bu_{t}$. The structural shocks $u_{t}\in R^{3}$ are
independent Gaussian with $\Lambda=diag(\lambda_{1},\lambda_{2},\lambda_{3})$
and $(\lambda_{1},\lambda_{2},\lambda_{3})=(2,0.4,0.2),(2,0.2,0.2),(1,1,1)$.
The structural change occurs at $\tau=0.4$ and the instantaneous effect matrix
$B$ is.

$B=\left[\begin{array}[]{ccc}27.92&0.231&1.569\\\ 0.441&5.643&0.079\\\
0.496&0.643&-4.668\\\ \end{array}\right].$

The choice of DGP4 is inspired by the empirical example used in the Section 5.

In order to enable a direct comparison of the MC results with the recent
outcomes presented by Lütkepohl et al. (2021), the same DGPs are selected for
analysis. This simulations are extended by DGP 4, which helps to evaluate
findings obtained during empirical application of the method.

In this novel testing procedure, the structural model is fitted to two sub-
samples of the size $T/2$, so the sample sizes used in MC are doubled as
compared to Lütkepohl et al. (2021). As the result, $T\in\\{200,500,1000\\}$
is considered. Finally, to examine the asymptotic properties of the test, the
results for $T=2000$ are also analyzed. The number of Monte Carlo replications
is set to equal 1000. The sample splitting is repeated for 100 times, when the
multiple split approach is used.

###  4.1 Non-Gaussian residuals

First, let us analyze the results of DGP1 for two different methods of
selecting the kurtosis value: $\kappa=3$ and $\kappa=\hat{\kappa}$. The fist
case corresponds to the assumption of residuals normality, whereas the second
one represents the case, in which the residuals are potentially non-normal.
Table 1 presents the outcomes for a single split (W) and a combined multiple
split (AveW) statistics. Columns 3-6 and 7-10 show the empirical size and the
empirical power of the test, respectively, for the nominal significance level
5%.

First, let us analyze the results for Gaussian shocks. It can be observed that
for short samples ($T=200$) the basic, $W$, test has a size exceeding 5%. As
the sample size increases, the empirical size converges to its nominal level.
Hence, the test is asymptotically valid. At the same time, the size of the
$AveW$ test is lower than for $W$ and falls to zero once the sample reaches
1000 observations. This reflects the fact that harmonic mean provides
conservative p-values (see Vovk and Wang (2020)). When the power of these
tests are considered, it can be noticed that the $W$ test has relatively poor
power (reaching 50% for $T=1000$). At the same time, the power $AveW$ is much
higher and exceeds 80% for $T=2000$.

The results are quite different, when the test is used for a non Gaussian
distribution. If the errors are $t$-distributed then their kurtosis differ
from $3$ and $W$ statistic is misspecified. As the result, $W$ is considerably
oversized even in large samples. For example, for $T=2000$, the test rejects
the null in more than 10% cases for 5% significance level. In contrast, when
$\kappa=\hat{\kappa}$ is allowed the rejection frequencies becomes much closer
to the significance level when the null hypothesis is true. For example, for
$T=1000$ the test reject $H_{0}$ in 7% cases and one could observe a slow
convergence to the nominal level. The power of the test is comparable with
that for Gaussian processes.

Finally, one can not observe large differences between Gaussian and t-Student
distributions for $AveW$ test. In case of non Gaussian errors, its size is
slightly larger than previously and exceeds 5%, but the difference diminish
once a longer sample ($T=500$) is used. Moreover, the power of the $AveW$
remains at a similar level even when the distribution changes. As the result,
the $AveW$ outperforms basic $W$ test in terms of the power level and the
robustness to the error distribution.

In the remaining part of the article, it is assumed that the kurtosis
parameter is estimated, $\kappa=\hat{\kappa}$.

Table 1: Empirical size and power of identification tests $W$ and $AveW$ ,
DGP1

| Size | Power
---|---|---
| $(\lambda_{1},\lambda_{2})=(2,1)$ | $(\lambda_{1},\lambda_{2})=(2,1)$
Kurtosis | Sample | Gaussian | $t(5)$ | Gaussian | $t(5)$
($\kappa$) | ($T$) | W | AveW | W | AveW | W | AveW | W | AveW
3 | 200 | 0.126 | 0.069 | 0.180 | 0.114 | 0.366 | 0.532 | 0.348 | 0.503
| 500 | 0.077 | 0.004 | 0.132 | 0.019 | 0.459 | 0.644 | 0.404 | 0.639
| 1000 | 0.045 | 0.000 | 0.100 | 0.001 | 0.487 | 0.745 | 0.490 | 0.723
| 2000 | 0.056 | 0.000 | 0.111 | 0.000 | 0.552 | 0.829 | 0.505 | 0.786
$\hat{\kappa}$ | 200 | 0.112 | 0.068 | 0.133 | 0.094 | 0.365 | 0.510 | 0.382 | 0.534
| 500 | 0.059 | 0.006 | 0.084 | 0.016 | 0.427 | 0.637 | 0.483 | 0.676
| 1000 | 0.055 | 0.000 | 0.072 | 0.000 | 0.494 | 0.740 | 0.532 | 0.760
| 2000 | 0.054 | 0.000 | 0.060 | 0.000 | 0.559 | 0.833 | 0.565 | 0.852

Note: for $AveW$: number of splits is $N=100$ and the truncation is at 20% and
80% percentiles

###  4.2 Estimation of VAR parameters

Next, let us analyze the impact of the estimation of VAR parameters on the
test performance. DGP2 is generated as a VAR(2) process, which parameters
correspond to the estimates of the Blanchard and Quah (1986) model fitted to
updated time series, as in Chen and Natšunajev (2016). As noticed by the
Lütkepohl et al. (2021), the performance of the identification test may depend
on the proper selection of the autoregression order. Therefore, the testing
procedure was conducted for both VAR(2) and VAR(1) specifications. The results
are presented in Table 2.

The outcomes of both $W$ and $AveW$ tests are similar to DGP1. One cannot
observe any significant change of the size or the power of these test when the
misspecified VAR(1) model is fitted. In case of both VAR(2) and VAR(1) models,
the size of the basic $W$ test is close to the nominal level, whereas the
$AveW$ test is too conservative. The power of the $W$ oscillates around 50%
for long samples. Finally, the power of $AveW$ test reaches 80% for $T=2000$.

Table 2: Empirical size and power of identification tests $W$ and $AveW$, DGP2

| Size | Power
---|---|---
| $(\lambda_{1},\lambda_{2})=(0.5,0.1)$ | $(\lambda_{1},\lambda_{2})=(0.5,0.5)$
Sample | VAR(1) | VAR(2) | VAR(1) | VAR(2)
($T$) | W | AveW | W | AveW | W | AveW | W | AveW
200 | 0.063 | 0.000 | 0.067 | 0.000 | 0.366 | 0.637 | 0.340 | 0.588
500 | 0.060 | 0.000 | 0.054 | 0.000 | 0.429 | 0.709 | 0.428 | 0.712
1000 | 0.036 | 0.000 | 0.052 | 0.000 | 0.494 | 0.788 | 0.484 | 0.758
2000 | 0.052 | 0.000 | 0.055 | 0.000 | 0.552 | 0.836 | 0.513 | 0.832

Note: kurtosis parameter is $\kappa=\hat{\kappa}$; for $AveW$: number of
splits is $N=100$ and the truncation is at 20% and 80% percentiles

###  4.3 Processes of higher dimension

In many empirical problems, the researchers need to analyze higher-dimension
processes. The approach proposed by Lütkepohl et al. (2021) experiences some
problems in case of three-dimensional Gaussian DGP due to a need of sequential
testing when a type of the identification violation is not known. The method
presented in this article solves some of the previous issues as it is
asymptotic valid regardless of the process dimension. Clearly, it is expected
that long samples are needed to estimate the model parameters accurately but
the asymptotic properties should be maintained under the null hypothesis.

The simulation results for the three dimensional case are presented in Table
3-4. In these tables, the results of $W$ and $AveW$ approaches are accompanied
by the outcomes of the $Q$ test proposed by Lütkepohl et al. (2021). It should
be notices that in case of $Q$, there are three potential null hypothesis
($\lambda_{1}=\lambda_{2}$, $\lambda_{2}=\lambda_{3}$,
$\lambda_{1}=\lambda_{2}=\lambda_{3}$). The $p$-value presented in the is
tables is the smallest $p$-value of the these three tests. Moreover, it should
be rememberd that the hypothesis of $Q$ test are stated differently: under the
null, the model is not identified whereas under the alternative it is. Hence,
the presented rejections of null have different interpretation: for $W$ and
$AveW$ tests they show how often the approach indicates lack of identification
and for $Q$ it shows the opposite.

For GDP3 the parameters $(\lambda_{1},\lambda_{2},\lambda_{3})=(3,2,1)$
represents the process under the null, whereas
$(\lambda_{1},\lambda_{2},\lambda_{3})=(2,2,2)$ or $(3,2,2)$ indicate the
model is not identifiable. For the DGP4, the null is described by
$(\lambda_{1},\lambda_{2},\lambda_{3})=(2,0.4,0.2)$. The results confirm that
the test is oversize for short samples. As expected, the empirical size
converges to the nominal level as the sample grows but substantially exceeds
5% for $T<2000$. When the two test $W$ and $AveW$ are considered, the results
supports previous findings. The multiple split approach is conservative as its
empirical size converges to zero. Moreover, it has substantially higher power
than the $W$ test, particularly for $(\lambda_{1},\lambda_{2},\lambda_{3})$
equal to $(3,2,2)$ for GDP3 and $(2,0.2,0.2)$ for GDP4.

Next, let us compare the presented approaches with the $Q$ statistic. It can
be noticed first that the Lütkepohl et al. (2021) test is too conservative as
its empirical size falls short the nominal 5%. Second, for very short samples,
$T=200$ it may have almost no power: it indicates lack of identification in
95.8% of cases when $(\lambda_{1},\lambda_{2},\lambda_{3})=(3,2,1)$. Finally,
as sample grows, the power of the test rises substantially and reaches almost
100% for the longest samples $T=2000$.

The results indicate that for short samples, both testing approaches are prone
to indicate falsely lack of identification. Hence, it should be used with some
concern. Despite this fact, the $AveW$ test can be useful when evaluating the
null even for $T=200$, as its size adjusted power is and 0.474 and 0.428 for
DGP3 and DGP4, respectively.

Table 3: Empirical size and power of tests $W$, $AveW$ and $Q$ for DGP3

Parameters | T | W | AveW | Q
---|---|---|---|---
(3,2,1) | 200 | 0.319 | 0.394 | 0.042
| 500 | 0.180 | 0.143 | 0.470
| 1000 | 0.113 | 0.038 | 0.824
| 2000 | 0.072 | 0.02 | 0.992
(3,2,2) | 200 | 0.543 | 0.868 | 0.000
| 500 | 0.548 | 0.852 | 0.000
| 1000 | 0.537 | 0.821 | 0.000
| 2000 | 0.567 | 0.829 | 0.000
(2,2,2) | 200 | 0.650 | 0.957 | 0.000
| 500 | 0.765 | 0.989 | 0.003
| 1000 | 0.843 | 0.994 | 0.013
| 2000 | 0.881 | 1.000 | 0.045

Note: kurtosis parameter is $\kappa=\hat{\kappa}$; for $AveW$: number of
splits is $N=100$ and the truncation is at 20% and 80% percentiles; for $Q$
test under $H_{0}$ the model parameters are not identified

Table 4: Empirical size and power of tests $W$, $AveW$ and $Q$ for DGP4

Parameters | T | W | AveW | Q
---|---|---|---|---
(2,0.4,0.2) | 200 | 0.148 | 0.087 | 0.560
| 500 | 0.071 | 0.003 | 0.938
| 1000 | 0.045 | 0.000 | 1.000
| 2000 | 0.054 | 0.000 | 1.000
(2,0.2, 0.2) | 200 | 0.369 | 0.515 | 0.048
| 500 | 0.422 | 0.623 | 0.050
| 1000 | 0.472 | 0.713 | 0.059
| 2000 | 0.536 | 0.813 | 0.039
(1,1,1) | 200 | 0.774 | 0.993 | 0.000
| 500 | 0.828 | 1.000 | 0.003
| 1000 | 0.846 | 0.998 | 0.013
| 2000 | 0.881 | 1.000 | 0.000

Note: kurtosis parameter is $\kappa=\hat{\kappa}$; for $AveW$: number of
splits is $N=100$ and the truncation is at 20% and 80% percentiles; for $Q$
test under $H_{0}$ the model parameters are not identified

###  4.4 Robustness analysis of $AveW$ approach

In the $AveW$ test, the $p$-values, $v_{n}$, stemming from different sample
splits are averaged via the harmonic mean. Before the mean is computed, the
extreme values of $v_{n}$ are removed. In the above results, 20% lowest and
highest observations are truncated. This rises a question, how the adopted
truncation threshold impacts the results. Therefore, different thresholds:
(10%,90%) and (30%, 70%) are applied to the DGP1. The results are presented in
Table 5, in which the first column shows the rejection frequency for the
specification previously studied. Next columns describe the outcomes of
different specifications. They indicate that the size and the power depends on
the truncation used. For (10%,90%) the null is rejected more often than for
(20%,80%) and (30%,70%). This leads to a growth of power at the cost of
exceeding the nominal size. When the adjusted power is analyzed, it shows that
widening of the truncation bounds may be preferable. Although the results are
quantitatively different, they reflect a similar pattern. Generally, the
$AveW$ test is conservative, its power increases with the sample size and is
larger than for $W$ test.

Another important decision, which needs to be taken when applying $AveW$ is
the choice of $N$: the number of $v_{n}$ to be averaged. In the above MC, $N$
is set to 100. In order to evaluate the robustness of the outcomes to the
selection of $N$, two other values are used: $N=50,200$. The resulting
rejections of the null hypothesis for DGP1 are presented in the last two
columns of Table 5. They show that the outcomes are robust. Change of $N$ and
does not alter significantly neither empirical size nor power.

To sum up, when applying the $AveW$ method, the researcher can choose any $N$,
as long as its big enough. It should however select the truncation thresholds
more carefully, taking into account that it may increase the size and result
in an over liberal approach, as for $T=200$ and (10%,90%).

Table 5: Rejection frequencies of the null hypothesis of $AveW$ for different
truncation thresholds and numbers the averaged $p$-values, $N$

| Truncation | (0.2,0.8) | (0.1,0.9) | (0.3,0.7) | (0.2,0.8) | (0.2,0.8)
---|---|---|---|---|---|---
$(\lambda_{1},\lambda_{2})$ | T/N | 100 | 100 | 100 | 50 | 200
(2,1) | 200 | 0.068 | 0.146 | 0.034 | 0.06 | 0.069
(size) | 500 | 0.006 | 0.018 | 0.000 | 0.004 | 0.002
| 1000 | 0.000 | 0.001 | 0.000 | 0.000 | 0.000
(2,2) | 200 | 0.532 | 0.694 | 0.377 | 0.497 | 0.522
(power) | 500 | 0.644 | 0.814 | 0.516 | 0.662 | 0.666
| 1000 | 0.745 | 0.877 | 0.597 | 0.732 | 0.747
(2,2) | 200 | 0.464 | 0.548 | 0.343 | 0.437 | 0.453
(adjusted power) | 500 | 0.638 | 0.796 | 0.516 | 0.658 | 0.664
| 1000 | 0.745 | 0.876 | 0.597 | 0.732 | 0.747

Note: kurtosis parameter is $\kappa=\hat{\kappa}$

##  5 Empirical example

The empirical example corresponds to the Kilian (2009) model of a global oil
market. In Kilian (2009) a VAR model is applied to describe the joint behavior
of three variables: $y_{t}=(\Delta prod_{t},q_{t},p_{t})$, where $\Delta
prod_{t}$ is the percentage growth of the global oil production, $q_{t}$ is
the logarithm of a detreded index of real economic activity and $p_{t}$ is the
logarithm of the real oil price. In the article, three structural shocks are
analyzed: oil supply shock, aggregated demand shock and an oil-market-specific
demand shock. In order to identify the shocks, the instantaneous effect
matrix, $B$, was assumed to be lower triangular.

The model has been evaluated by Lütkepohl and Natšunajev (2014), who explore
heteroscedasticity to identify the structural parameters. In the paper, the
monthly data from 1973m2–2006m12 are used to estimate a Markow Switching
model. The results indicate that there are only weak evidence of a lower
triangular form of the matrix $B$. This conclusion was next challenged by
Lütkepohl et al. (2021), where the shift of shocks volatility in October 1987
was used to identify the model. The outcomes suggest that the Kilian model can
be only partially identified via the variance change and hence the lower-
triangularity test should be compared with $\chi^{2}(2)$ distribution, rather
than $\chi^{2}(3)$. These results show that identification testing can be
crucial for interpreting the model outcomes.

Here we use the data explored by Lütkepohl et al. (2021) and apply the
proposed testing procedure to verify if the change of shocks volatility is
sufficient to indentify the model structure. Similar to previous papers, a
VAR(3) model is fitted to the data with the GLS method. The estimated relative
variances are $\lambda_{1}=0.186$, $\lambda_{2}=0.363$ and
$\lambda_{3}=2.157$. The full samples and sub-samples estimator of the
randomly split sample are presented in Table 6. Although the estimators
$\hat{\theta}_{1}$ and $\hat{\theta}_{2}$ are distinct, they are very similar,
which supports the null hypothesis.

Table 6: Estimates of model parameters for the whole sample ($\hat{\theta}$) and a randomly split sub-samples ($\hat{\theta}_{1}$, $\hat{\theta}_{2}$). Parameter | $\hat{\theta}$ | $\hat{\theta}_{1}$ | $\hat{\theta}_{2}$
---|---|---|---
$B_{11}$ | 27.918 | 30.493 | 24.879
$B_{12}$ | 0.441 | 0.233 | 0.539
$B_{13}$ | 0.496 | 0.267 | 0.694
$B_{21}$ | 0.231 | -0.825 | 1.756
$B_{22}$ | 5.643 | 5.258 | 6.003
$B_{23}$ | 0.643 | 0.089 | 1.137
$B_{31}$ | 1.569 | 1.203 | 1.792
$B_{32}$ | 0.079 | 0.028 | 0.122
$B_{33}$ | -4.668 | -3.378 | -5.624
$\lambda_{1}$ | 0.186 | 0.184 | 0.189
$\lambda_{2}$ | 0.363 | 0.404 | 0.332
$\lambda_{3}$ | 2.157 | 4.173 | 1.461

Finally, the $W_{n}$ statistics are computed and compared with $\chi^{2}(12)$
distribution. The individual and averaged p-values are $v_{1}=0.5096$ and
$\bar{v}=1$, which indicate that there are no strong evidences against the
null hypothesis. The results support initial conclusions of Lütkepohl and
Natšunajev (2014) and question the lower-traingularity of the $B$ matrix.
However, it should be noticed here that the Monte Carlo results show that the
size adjusted power of the test for $T=500$ is 0.620, which is much below 1.
Therefore, the conclusions presented are not very strong and should be
repeated once more data becomes available.

##  6 Conclusions

In this article, a novel test is proposed to verify if the parametric model is
identifiable. Unlike in previous papers, the test null hypothesis assumes
model identification, whereas under the alternative the model parameters are
not identified. This modifications allows to use the well known results of the
ML estimation method, which states that under some regularity conditions ML
estimators are consistent and asymptotically normal. In order to explore these
features, two parameter estimators are compared, based on a random split of
the sample. Under the null, both estimators are consistent and hence the
distance between them converges to zero. The significance of the differences
is assessed with a Wald type of test, which has a well known $\chi^{2}$
distributions with the number of degrees of freedom corresponding to the
number of parameters. Next a modification of the approach is proposed, which
explores outcomes from multiple splitting of the sample. The resulting
$p$-values are combined with a truncated, scaled harmonic mean.

Next, the method is adjusted for testing whether heteroscedasticity is
sufficient to identify parameters of a SVAR model. In order to assess the
properties of the test, a MC experiment is conducted. The data generating
processes used here correspond to those applied by Lütkepohl et al. (2021).
The results of the analysis indicate that for short samples, the $W$ test is
slightly over-sized. However, as the sample size increases, the empirical size
converges to its nominal level. The power of the test depends on the dimension
of the VAR model. It improves as the number of identical variances increases.
It is also demonstrated that the multiple split test, $AveW$, is conservative,
with its empirical size falling short of the nominal level. At the same time,
the power to $AveW$ exceeds the power of $W$ counterpart. Moreover, it is more
robust to VAR model misspecification and non-Gaussian distribution of
residuals. These results shows that $AveW$ is a reasonable choice in empirical
applications.

I believe that the theory can be developed in various directions. First, it
can be applied to SVAR models with more complex heteroscedasticity structure,
as in Maciejowska et al. (2016), Lütkepohl and Natšunajev (2014). In such case
a different sampling method should be considered. Moreover, it can be used to
assess many other models, such as MN, MS or STR, which are affected by the
identification problem. It would be desirable to evaluate the empirical
properties of the test in the context of a broad family of applications.

## References

  * Andrews (1993) Andrews, D.W.K., 1993.  Tests for parameter instability and structural change with unknown change point.  Econometrica 61, 821–856.
  * Andrews and Ploberger (1994) Andrews, D.W.K., Ploberger, W., 1994\.  Optimal tests when a nuisance parameter is present only under the alternative.  Econometrica 62, 1383–1414.
  * Blanchard and Quah (1986) Blanchard, Q.J., Quah, D., 1986\.  The dynamic effects of aggregate demand and supply disturbances.  American Economic Review 34, 523–533.
  * Chen and Natšunajev (2016) Chen, W., Natšunajev, A., 2016\.  On the long-run neutrality of demand shocks.  Economics Letters 139, 57–60.
  * van Dijk et al. (2002) van Dijk, D., Teräsvirta, T., Franses, P.H., 2002.  Smooth Transition Autoregressive models - a survey of recent developments.  Econometric Reviews 21, 1–47.
  * Greene (2018) Greene, W.H., 2018.  Econometric Analysis, 8th Edition.  Pearson Education, New Jersey.
  * Hansen (1992) Hansen, B.E., 1992.  The Likelihood Ratio test under nonstandard conditions: testing the Markow Switching model of GNP.  Journal of Applied Econometrics 7, 61–82.
  * Hansen (1996) Hansen, B.E., 1996.  Inference when a nuisance parameter is not identified under the null hypothesis.  Econometrica 64, 413–430.
  * Hsiao (1983) Hsiao, C., 1983.  Identification.  Griliches, Z. and M.D. Intriligator, (Eds.), Handbook of econometrics, Vol. 1. North Holland, Amsterdam, The Netherlands.
  * Kikian and Lütkepohl (2017) Kikian, L., Lütkepohl, H., 2017\.  Structural Vector Autoregressive Analysis.  Cambridge, UK: Cambridge University Press.
  * Kilian (2009) Kilian, L., 2009.  Not all oil shocks are alike: disentangling demand and supply shocks in the crude oil market.  American Economic Review 99, 1053–1069.
  * Lanne and Lütkepohl (2008) Lanne, M., Lütkepohl, H., 2008\.  Identifying monetary policy shocks via changes in volatility.  Journal of Money, Credit and Banking 40, 1131–1149.
  * Lanne et al. (2010) Lanne, M., Maciejowska, K., Lütkepohl, H., 2010.  Structural vector autoregressions with markov switching.  Journal of Economic Dynamics and Control 34, 121–131.
  * Lütkepohl (2005) Lütkepohl, H., 2005.  New Introduction to Multiple Time Series Analysis.  Springer-Verlag, Berlin.
  * Lütkepohl et al. (2021) Lütkepohl, H., Meitz, M., Netsunajev, A., Saikkonen, P., 2021\.  Testing identification via heteroskedasticity in structural vector autoregressive models.  Econometrics Journal 24, 1–214.
  * Lütkepohl and Natšunajev (2014) Lütkepohl, H., Natšunajev, A., 2014\.  Disentangling demand and supply shocks in the crude oil market: how to check sign restrictions in structural VARs.  Journal of Applied Econometrics 29, 479–496.
  * Lütkepohl and Natšunajev (2017) Lütkepohl, H., Natšunajev, A., 2017\.  Structural vector autoregressions with smooth transition in variances. The interaction between U.S. monetary policy and the stock market.  Econometrics and Statistics 1, 2–18.
  * Lütkepohl and Velinov (2021) Lütkepohl, H., Velinov, A., 2021\.  Structural vector autoregressions: checking identifying long-run restrictions via heteroscedasticity.  Journal of Economic Surveys 30, 3777–392.
  * Lütkepohl and Woźniak (2021) Lütkepohl, H., Woźniak, T., 2021\.  Bayesian inference for structural vector autoregressions identified by Markov-switching heteroskedasticity.  Journal of Economic Dynamics and Control 133, 103862.
  * Luukkonen et al. (1988) Luukkonen, R., Saikkonen, P., Ter asvirta, T., 1988.  Testing linearity against smooth transition autoregressive models.  Biometrika 75, 491–499.
  * Maciejowska et al. (2016) Maciejowska, K., Nowotarski, J., Weron, R., 2016.  Probabilistic forecasting of electricity spot prices using factor quantile regression averaging.  International Journal of Forecasting 32, 957–965.
  * McLachlan et al. (2018) McLachlan, G.J., Lee, S.X., Rathnayake, S.I., 2018.  Finite mixture models.  Annual Review of Statistics and Its Application 6, 355–378.
  * Rossi and Inoue (2011) Rossi, B., Inoue, A., 2011.  Testing for weak identification in possibly nonlinear models.  Journal of Econometrics 161, 246–261.
  * Rothenberg (1971) Rothenberg, T., 1971.  Identification in parametric models.  Econometrica 39, 577–591.
  * Schott (2001) Schott, J.R., 2001.  Some tests for the equality of covariance matrices.  Journal of Statistical Planing and Inference 94, 25–46.
  * Sims (1980) Sims, C.A., 1980.  Macroeconomics and reality.  Econometrica 48, 1–48.
  * Vovk and Wang (2020) Vovk, V., Wang, R., 2020.  Combining p-values via averaging.  Biometrica 7107, 1791–808.