# Gravitational waves from fragmentation of a primordial scalar condensate
into Q-balls

Alexander Kusenko  Department of Physics and Astronomy, University of
California, Los Angeles, CA 90095-1547, USA  Anupam Mazumdar  Physics
Department, Lancaster University, Lancaster, LA1 4YB, UK  Niels Bohr
Institute, Blegdamsvej-17, Copenhagen, DK-2100, Denmark

###### Abstract

A generic consequence of supersymmetry is formation of a scalar condensate
along the flat directions of the potential at the end of cosmological
inflation. This condensate is usually unstable, and it can fragment into non-
topological solitons, Q-balls. The gravitational waves produced by the
fragmentation can be detected by Laser Interferometer Space Antenna (LISA),
Advanced Laser Interferometer Gravitational-Wave Observatory (LIGO), and Big
Bang Observer (BBO), which can offer an important window on the early universe
and the physics at some very high energy scales.

###### pacs:

04.30.-w,04.30.Tv,11.30.Pb,12.60.Jv,98.80.-k

††preprint: UCLA/08/TEP/25

Supersymmetry is widely regarded as a likely candidate for physics beyond the
Standard Model. While many variants of supersymmetry have been considered, all
of them have scalar potentials with some flat directions lifted only by the
supersymmetry-breaking terms. At the end of cosmological inflation, the
formation of a scalar condensate along the flat directions can have a number
of important consequences reviews . In particular, it can be responsible for
generation of the matter-antimatter asymmetry via Affleck–Dine (AD) mechanism
AD , and, in some models, dark matter can be produced in the same process KS ;
dark_matter . Some flat directions could be responsible for the primordial
inflation Randall:1995dj ; AEGM .

The formation of AD condensate is a generic phenomenon, relying only on the
assumptions of inflation and supersymmetry. In general, this condensate is
unstable: an initially homogeneous condensate can break up into lumps of the
scalar field, called Q-balls Coleman , under some very generic conditions KS .
All phenomenologically acceptable supersymmetric generalizations of the
Standard Model admit Q-balls Kusenko:1997zq , which can be stable, or can
decay into fermions Kusenko:1997zq ; DKS . The formation of Q-balls is
accompanied by a coherent motion of the scalar condensate, which creates the
source of gravity waves. We will show that fragmentation of the scalar
condensate into Q-balls can produce gravitational waves detectable by LISA
LISA , LIGO III LIGOIII , and BBO BBO .

The physics of AD condensate fragmentation has been studied both analytically
KS ; EMc ; EJM ; Pawl ; Johnson:2008se  and numerically Minos ; KK ; MV ;
Volkov ; Palti ; Broadhead ; Campanelli . At the end of inflation (assuming
that the inflation occurs in a hidden sector at high scales) the condensate
has a uniform density, with small perturbations of the order of $10^{-5}$
reviews . Under some rather generic conditions, the instabilities develop and
lead, eventually, to the formation of Q-balls, which can either decay or
remain as stable relics KS ; EMc . Although the final product of such
evolution, Q-balls in the ground state, are spherically symmetric, the
coherent motions associated with the condensate fragmentation and re-
arrangement are not spherically symmetric. Moreover, the newly formed Q-balls
first appear in their excited states and oscillate until they settle in the
spherically symmetric ground states KK ; MV . The lack of spherical symmetry
in the process of fragmentation is essential for generating the gravity waves.

Following the general picture developed in Refs. KS ; EJM , the scalar
condensate undergoing fragmentation can be approximated, in the linear regime,
as $\phi(x,t)=\phi(t)\equiv R(t)e^{i\Omega(t)}$, plus a perturbation $\delta
R,\delta\Omega\propto e^{S(t)-i\vec{k}\vec{x}}$. One finds that the
homogeneous solution is unstable due to some exponentially growing modes,
${\rm Re}\,\alpha>0$, where $\alpha=dS/dt$ KS ; EJM ; Johnson:2008se . The
mass density of the condensate undergoing fragmentation can be written as

$\rho(x,t)=\rho_{0}+\rho_{1}(x,t)\,,$ (1)

where

$\rho_{1}(x,t)=\epsilon\rho_{0}\int d^{3}k\,e^{\alpha_{k}t}\cos(\omega
t-\vec{k}\cdot\vec{x})\,.$ (2)

The instability develops when there is a band of growing modes with positive
and large enough $\alpha_{k}$ KS . The linear approximation breaks down when
$\epsilon\exp(\alpha_{k}t)\sim 1$, but we will use this representation, up to
its limit of applicability, to get the estimates of the gravity waves
produced.

The quadrupole moment that generates gravity waves is given by Weinberg

$D_{ij}=\int d^{3}x\ x_{i}x_{j}\,T^{00}(x,t)\,,$ (3)

where the energy-momentum tensor $T^{00}(x,t)\approx\rho(x,t)$. The space
integration is over some arbitrary volume.

The power emitted in gravity waves in one frequency mode is given by:

$P(\omega)=\frac{2}{5}G\omega^{6}\left(D^{*}_{ij}(\omega)D_{ij}(\omega)-\frac{1}{3}|D_{ij}(\omega)|^{2}\right)\,,$
(4)

and the total energy emitted in gravitational waves, in all frequencies, is
given by:

$E\sim\left(2\pi\int p(\omega)d\omega\right)\times\Delta t\,,$ (5)

where $\Delta t$ is the duration of the fragmentation.

Based on the analytical and numerical calculations of the condensate
fragmentation KS ; KK ; Kasuya2 , we take the typical parameters of the
fastest-growing mode:

$k\sim\xi_{k}\times 10^{2}H_{\ast},\ \omega_{k}\sim vk\sim\xi_{k}\times
10^{2}\,vH_{\ast},$ (6)

where $H_{\ast}$ is the Hubble constant at the time of the condensate
fragmentation, and $v$ is the typical group velocity of the wave front in the
evolution of the condensate, and we expect that the dimensionless factor
$\xi_{k}\sim 1$, based on the results of Refs. KS ; KK ; Kasuya2 .

Since no cancellations are expected in the absence of spherical symmetry, we
replace the $x_{i}x_{j}$ by $(f_{k}\times 10^{2}H_{\ast})^{-2}$ in the space
integration, take the volume to be $V\sim H^{-3}_{\ast}$, and assume that
$\epsilon\exp(\alpha_{k}t)\sim 1$. Then, for the leading mode,

$D_{ij}(t)\sim
H_{\ast}^{-3}\left(10^{2}H_{\ast}\right)^{-2}\rho_{0}\cos(\omega_{k}t-kx)\,,$
(7)

and, in frequency space,

$D_{ij}(\omega)\sim 10^{-4}\xi_{k}^{-2}\frac{\rho_{0}}{H_{\ast}^{5}}\,.$ (8)

For $\omega\sim 10^{2}vH_{\ast}$, we estimate the power in gravitational waves
in a Hubble volume:

$P\sim 10^{4}\xi_{k}^{-2}\,G\frac{\rho_{0}^{2}v^{6}}{H^{4}_{\ast}}\,.$ (9)

To estimate the velocity of the wavefront in the process of fragmentation, we
note that, for the mode $\phi(x,t)\approx
R(t)\exp\\{\alpha_{k}t\\}\cos(\omega_{k}t-kx)$, where $R(t)$ is a slowly
changing function of time,111The adiabatic limit $\dot{R}/R\rightarrow 0$ is
amenable to perturbation theory KS . This case corresponds to a large global
charge density. For a rapidly varying $R(t)$, the numerical calculations give
similar results regarding the fragmentation time and length scales KK .

$v\sim\xi_{v}\left|\dot{\phi}/\phi^{\prime}_{x}\right|\sim\xi_{v}\
\alpha_{k}/k,\ \ \xi_{v}\sim 1,$ (10)

where the uncertainty factor $\xi_{v}\sim 1$ will be retained to keep track of
the uncertainty in the final answer.

The relation between $\alpha_{k}$ and $k$ is given by a dispersion relation KS
, which takes a simple form,

$(\alpha_{k}^{2}+k^{2})\left(\alpha_{k}^{2}-\left(\dot{\Omega}^{2}-V^{\prime\prime}(R)\right)\right)+4\dot{\Omega}^{2}\alpha_{k}^{2}=0,$
(11)

under the following assumptions: $H\ll
k\sim\alpha_{k}\ll\sqrt{\dot{\Omega}^{2}-V^{\prime\prime}(R)}\sim\dot{\Omega}\sim
m_{\phi}$ valid in the case of the fastest-growing mode in gravity mediated
supersymmetry breaking models (the latter is essential for the
$\dot{\Omega}\sim m_{\phi}$ condition reviews ). This equation has an
approximate solution: $\alpha_{k}\approx k/\sqrt{3}$, and

$v^{6}\sim\xi_{v}^{6}(\alpha_{k}/k)^{6}\sim\xi_{v}^{6}(1/\sqrt{3})^{6}\sim
10^{-2}\xi_{v}^{6}\,.$ (12)

The fragmentation takes place on the time scale of the order of $\Delta
t\sim\alpha_{k}^{-1}\sim\xi_{k}^{-1}10^{-2}H^{-1}_{\ast}$. (Here we neglect
the possible contributions from the collisions and oscillations of Q-balls,
which can take place on a much longer time scale KS ; MV1 .) The total energy
in gravity waves generated in the Hubble volume is

$E\sim P\Delta t\sim
G\frac{\rho_{0}^{2}}{H_{\ast}^{5}}\,\xi_{k}^{-3}\xi_{v}^{6}\,.$ (13)

This corresponds to the energy density in gravitational waves at the time of
production

$\rho_{GW\ast}\sim
10^{-3}\xi_{k}^{-3}\xi_{v}^{6}\,\frac{\rho_{0}^{2}}{H^{2}_{\ast}M_{\rm
Pl}^{2}}\,,$ (14)

where $M_{\rm Pl}=1/\sqrt{8\pi G}$ is the reduced Planck mass. Hence, the
fraction of the energy density in gravitational waves at the time of
production is

$\Omega_{GW*}\sim
10^{-3}\xi_{k}^{-3}\xi_{v}^{6}\,\frac{\rho_{0}^{2}}{(H_{\ast}M_{\rm
Pl})^{4}}\,.$ (15)

If the energy density of the condensate is comparable to the total energy
density, or if the condensate energy dominates the energy in the universe,
then $\rho_{0}\sim 3H_{\ast}^{2}M_{\rm Pl}^{2}$, and $\Omega_{GW*}\sim
10^{-3}$.

The energy density in the condensate depends on the model, and, foremost, on
the type of supersymmetry breaking terms that lift the flat direction. This is
because the potential along the flat direction depends on supersymmetry
breaking (it vanishes in the limit of exact supersymmetry), and there are many
ways to break supersymmetry. In gauge-mediated supersymmetry breaking
scenarios the potential can have the form KS

$V(\phi)\approx M_{S}^{4}\log\left(1+\frac{|\phi|^{2}}{M_{S}^{2}}\right)\,.$
(16)

Here $M_{S}$ is the scale of supersymmetry breaking, which is of the order of
${\cal O}(1)$ TeV. In gravity mediated scenarios, the flat directions are
lifted by mass terms that persist all the way to the Planck scale EMc :

$V(\varphi)\approx
m_{\phi}^{2}\left(1+K\log\left(\frac{|\phi|^{2}}{M_{Pl}^{2}}\right)\right)|\phi|^{2}\,,$
(17)

where $K\sim 0.05$ (for squark directions) describes the running of the mass
term EMc . Since $m_{\phi}\sim M_{S}\sim 1-10$ TeV in typical models, both
potentials are phenomenologically acceptable near the minimum. However, in AD
condensate and inside the Q-balls that form in its fragmentation, the vacuum
expectation value (VEV) can be very large.

The main difference between gauge and gravity mediated cases for us is the
mass per baryon number stored in the AD condensate and in the Q-balls that
form eventually as a result of the fragmentation. In the gravity mediated
scenarios, the mass density is $\rho_{0}\sim m_{\phi}^{2}\phi^{2}$, the global
charge density is $n_{Q}\sim m_{\phi}\phi^{2}$, and the mass per unit global
charge is of the order of $m_{\phi}$, independent of the VEV $\phi_{0}$. In
gauge mediated scenarios, the mass density is $\rho_{0}\sim m_{\phi}^{4}$, the
global charge density is $n_{Q}\sim m_{\phi}\phi^{2}$, and the mass per unit
global charge is $\rho_{0}/n_{Q}\sim m_{\phi}^{2}/\phi$ KS ; DKS .

The flat directions that carry a non-zero global charge $Q=(B-L)$ contribute
to the generation of the baryon asymmetry via AD process AD . The requirement
that $\eta_{B}=n_{B}/n_{\gamma}\sim 10^{-10}$ implies that the total mass
density of such a condensate cannot be of the order of the total density of
the universe in generic models. This is true in both gauge-mediated and
gravity-mediated supersymmetry breaking models. The gravity waves from the
fragmentation of such a condensate are well below the capabilities of the
current and planned detectors.

On the other hand, there are flat directions whose baryon number $B$ and
lepton number $L$ are equal to each other DRT ; GKM ; reviews . While the flat
directions with $B\neq L$ contribute to baryon asymmetry of the universe,
those with $B=L$, have zero $(B-L)$ density. Electroweak sphalerons destroy
any primordial $(B+L)$ asymmetry, and so the corresponding
$\eta_{B+L}=n_{B+L}/n_{\gamma}$ is not constrained. It is possible that, at
the time of the fragmentation, $\eta_{B+L}\gg\eta_{B}$. For $B=L$ flat
directions, there is no reason why $\rho_{0}$ cannot be of the order of the
total energy density. The fragmentation of such flat directions can produce a
detectable level of gravitational waves.

There are various examples in the literature of the flat directions that can
dominate the energy density of the universe, while they do not contribute to
(and are not constrained by) the observed baryon asymmetry of the universe.
This is the case, for example, when the effective mass for the phase direction
is large during inflation, which results in the initial condition with a very
small $\dot{\Omega}$ for the scalar condensate $\phi(x,t)=\phi(t)\equiv
R(t)e^{i\Omega(t)}$. This is also the case when the inflation is driven by a
flat direction, $udd$, $LLe$ or $NH_{u}L$ Randall:1995dj ; AEGM ; AEGJM .
Another well-know example is a flat direction that acts as a curvaton and
dominates the energy density of the universe at the time of oscillations and
decay Curvaton . In all of these cases, the net global charge of the
condensate is negligible, but the fragmentation can still occur and produce
Q-balls and anti-Q-balls Kasuya1 ; Kasuya2 .

Once the gravitational waves are created, they are decoupled from the rest of
the plasma. We can estimate the peak frequency of the gravitational radiation
observed today, $f_{\ast}=\omega_{k}/2\pi$:

$\displaystyle f$ $\displaystyle=$ $\displaystyle
f_{\ast}\frac{a_{\ast}}{a_{0}}=f_{\ast}\left(\frac{a_{\ast}}{a_{\rm
rh}}\right)\left(\frac{g_{s,0}}{g_{s,{\rm
rh}}}\right)^{1/3}\left(\frac{T_{0}}{T_{\rm rh}}\right)\,$
$\displaystyle\approx$ $\displaystyle 0.6\,{\rm mHz}\
\xi_{k}\xi_{v}\,\left(\frac{g_{s,{\rm
rh}}}{100}\right)^{1/6}\left(\frac{T_{\rm rh}}{1~{}{\rm
TeV}}\right)\left(\frac{f_{\ast}}{10H_{\ast}}\right)\,,$

where we have assumed that $a_{\ast}\approx a_{\rm rh}$, which also means that
during the oscillations of the AD condensate the effect of the Hubble
expansion is negligible. The values of relativistic degrees of freedom are
$g_{s,{\rm rh}}\sim 300$, $g_{s,0}\sim 3.36$. The subscript “rh” denotes the
epoch of reheating and thermalization, while the subscript “0” refers to the
present time. As we discussed, the typical frequency of the oscillations of
the AD condensate is $\omega_{k}\sim 10^{2}H_{\ast}$ KS ; KK ; Kasuya2 . Then
for $T_{\rm rh}\sim 1$ TeV (such a value of the reheat temperature is natural
when the flat direction is responsible for reheating the universe Rouzbeh ;
AEGJM ), the frequency is of the order of mHz, which is in the right frequency
range for LISA LISA . A higher temperature $T_{\rm rh}\sim 100$ TeV
corresponds to the LIGO III frequency range, $10-100$ Hz LIGOIII . Signals in
both of these ranges will be accessible to BBO BBO . Since the supersymmetry
breaking scale is related to the energy in the condensate, as well as the
reheating temperature, LIGOIII and BBO could be in the position to probe
supersymmetry broken above 100 TeV, beyond the reach of Large Hadron Collider
(LHC).

The fraction of the critical energy density $\rho_{c}$ stored in the gravity
waves today is

$\displaystyle\Omega_{GW}$ $\displaystyle=$
$\displaystyle\Omega_{GW^{\ast}}\left(\frac{a_{\ast}}{a_{0}}\right)^{4}\left(\frac{H_{\ast}}{H_{0}}\right)^{2}$
$\displaystyle\approx$ $\displaystyle\frac{1.67\times
10^{-5}}{h^{2}}\left(\frac{100}{g_{s,\ast}}\right)^{1/3}\Omega_{GW^{\ast}}\approx
10^{-8}\,\xi_{k}^{-3}\xi_{v}^{6}\,h^{-2}$

where $a_{0}$ and $H_{0}$ are the present values of the scale factor and the
Hubble expansion rate. LISA can detect the gravitational waves down to
$\Omega_{GW}h^{2}\sim 10^{-11}$ at mHz frequencies, while LIGO III is
sensitive to $\Omega_{GW}h^{2}\sim(10^{-5}-10^{-11})$ in the $(5-10^{3})$ Hz
frequency band. Therefore, the gravitational waves with $\Omega_{GW*}\sim
10^{-3}$ and $\Omega_{GW}h^{2}\sim 10^{-8}$ from the fragmentation of an AD
condensate can be detected. The first results of our numerical simulations
(work in progress) appear to produce the gravitational wave signal that is
somewhat weaker. We attribute the difference to the value of the uncertainty
factor $\xi_{k}^{-3}\xi_{v}^{6}$, which is especially sensitive to the average
wavefront velocity.

As one can see from eq. (14), the power generated in frequency $\omega$ is
proportional to

$\xi_{k}^{-3}\sim\left(\frac{\omega}{10^{2}H_{*}}\right)^{-3}.$

Hence, the spectrum is strongly peaked near the longest wavelength, of the
order of the Q-ball size KS ; KK ; Kasuya2 . The relatively narrow spectral
width will help distinguish this signal from the gravity waves generated by
inflation inflation , which are expected to have an approximately scale-
invariant spectrum (and a smaller amplitude). Future numerical simulations
will help refine the prediction for the signal from Q-ball formation, which
can help distinguish this source from a phase transition in the early universe
Kamionkowski:1993fg . LISA and LIGO III will be able to discriminate the
gravity waves due to fragmentation from those of point sources, such as
merging black holes and neutron stars, which have specific “chirp” properties
Owen:1998dk . Furthermore, the signal discussed here will not create a
significant background for the cosmic microwave polarization experiments, such
as B-Pol B-pol , which can detect the gravity waves with extremely long
wavelength.

Some additional gravitational waves can be generated by collisions and
oscillations of Q-balls KS ; MV1 . We leave the discussion of the magnitude of
this additional contribution to future studies.

To summarize, the fragmentation of a scalar condensate into Q-balls, which is
a generic consequence of supersymmetry and inflation, can produce a detectable
level of gravitational waves, up to $\Omega_{GW}h^{2}\sim 10^{-8}$, near the
peak frequency of BOO and either LISA or LIGO, depending on the reheating
temperature. Detection of the gravitational waves form this process can shed
light on the earliest post-inflationary epoch in the history of the universe,
can probe supersymmetry even if it is broken at a scale above 100 TeV, and can
provide information about new physics at some very high energy scales
associated with the flat directions.

A.K. thanks S. Phinney for helpful discussions. The work of A.K. was supported
in part by the DOE grant DE-FG03-91ER40662 and by the NASA ATFP grant
NNX08AL48G. Research of A.M. is supported in part by the grant MRTN-
CT-2006-035863. A.K. thanks Aspen Center for Physics for hospitality.

## References

  * (1) K. Enqvist and A. Mazumdar, Phys. Rept. 380, 99 (2003); M. Dine and A. Kusenko, Rev. Mod. Phys. 76, 1 (2004).
  * (2) I. Affleck and M. Dine, Nucl. Phys. B 249, 361 (1985).
  * (3) A. Kusenko and M. E. Shaposhnikov, Phys. Lett. B 418, 46 (1998).
  * (4) A. Kusenko, V. Kuzmin, M. E. Shaposhnikov and P. G. Tinyakov, Phys. Rev. Lett. 80, 3185 (1998); M. Laine and M. E. Shaposhnikov, Nucl. Phys. B 532, 376 (1998); A. Kusenko, L. Loveridge and M. Shaposhnikov, Phys. Rev. D 72, 025015 (2005); JCAP 0508, 011 (2005).
  * (5) L. Randall, M. Soljacic and A. H. Guth, Nucl. Phys. B 472, 377 (1996).
  * (6) R. Allahverdi, K. Enqvist, J. Garcia-Bellido and A. Mazumdar, Phys. Rev. Lett. 97, 191304 (2006). R. Allahverdi, A. Kusenko and A. Mazumdar, JCAP 0707, 018 (2007).
  * (7) G. Rosen, J. Math. Phys. 9 (1968) 996; 9 (1968) 999; R. Friedberg, T. D. Lee and A. Sirlin, Phys. Rev. D 13, 2739 (1976); S. R. Coleman, Nucl. Phys. B 262, 263 (1985) [Erratum-ibid. B 269, 744 (1986)]; A. Kusenko, Phys. Lett. B 404, 285 (1997).
  * (8) A. Kusenko, Phys. Lett. B 405, 108 (1997);
  * (9) G. R. Dvali, A. Kusenko and M. E. Shaposhnikov, Phys. Lett. B 417, 99 (1998).
  * (10) http://lisa.nasa.gov/
  * (11) http://www.ligo.caltech.edu/
  * (12) G. M. Harry, et al., Class. Quant. Grav. 23, 4887 (2006) [Erratum-ibid. 23, 7361 (2006)]; V. Corbin and N. J. Cornish, Class. Quant. Grav. 23, 2435 (2006); http://universe.nasa.gov/new/program/bbo.html
  * (13) K. Enqvist and J. McDonald, Phys. Lett. B 425, 309 (1998). K. Enqvist and J. McDonald, Nucl. Phys. B 538, 321 (1999).
  * (14) A. Pawl, Nucl. Phys. B 679, 231 (2004).
  * (15) K. Enqvist, A. Jokinen and J. McDonald, Phys. Lett. B 483, 191 (2000).
  * (16) M. C. Johnson and M. Kamionkowski, arXiv:0805.1748 [astro-ph].
  * (17) M. Axenides, S. Komineas, L. Perivolaropoulos and M. Floratos, Phys. Rev. D 61, 085006 (2000).
  * (18) S. Kasuya and M. Kawasaki, Phys. Rev. D 61, 041301 (2000); S. Kasuya and M. Kawasaki, Phys. Rev. D 62, 023512 (2000). S. Kasuya and M. Kawasaki, Phys. Rev. Lett. 85, 2677 (2000). S. Kasuya and M. Kawasaki, Phys. Rev. D 64, 123515 (2001).
  * (19) T. Multamaki and I. Vilja, Nucl. Phys. B 574, 130 (2000). K. Enqvist, A. Jokinen, T. Multamaki and I. Vilja, Phys. Rev. D 63, 083501 (2001). T. Multamaki and I. Vilja, Phys. Lett. B 535, 170 (2002). T. Multamaki, Phys. Lett. B 511, 92 (2001).
  * (20) M. S. Volkov and E. Wohnert, Phys. Rev. D 66, 085003 (2002).
  * (21) E. Palti, P. M. Saffin and E. J. Copeland, Phys. Rev. D 70, 083520 (2004). M. I. Tsumagari, E. J. Copeland and P. M. Saffin, arXiv:0805.3233 [hep-th].
  * (22) M. Broadhead and J. McDonald, Phys. Rev. D 72, 043519 (2005). M. Broadhead and J. McDonald, Phys. Rev. D 69, 063510 (2004).
  * (23) L. Campanelli and M. Ruggieri, Phys. Rev. D 77, 043504 (2008).
  * (24) S. Weinberg, Gravitation and Cosmology, John Wiley $\&$ Sons. 1972.
  * (25) K. Enqvist, S. Kasuya and A. Mazumdar, Phys. Rev. D 66, 043505 (2002).
  * (26) T. Multamaki and I. Vilja, Phys. Lett. B 482, 161 (2000). T. Multamaki and I. Vilja, Phys. Lett. B 484, 283 (2000).
  * (27) K. Enqvist, S. Kasuya and A. Mazumdar, Phys. Rev. Lett. 89, 091301 (2002).
  * (28) M. Dine, L. Randall and S. D. Thomas, Nucl. Phys. B 458, 291 (1996).
  * (29) T. Gherghetta, C. F. Kolda and S. P. Martin, Nucl. Phys. B 468, 37 (1996).
  * (30) R. Allahverdi, K. Enqvist, J. Garcia-Bellido, A. Jokinen and A. Mazumdar, JCAP 0706, 019 (2007).
  * (31) K. Enqvist, S. Kasuya and A. Mazumdar, Phys. Rev. Lett. 90, 091302 (2003). M. Postma, Phys. Rev. D 67, 063518 (2003). K. Enqvist, A. Jokinen, S. Kasuya and A. Mazumdar, Phys. Rev. D 68, 103507 (2003). R. Allahverdi, K. Enqvist, A. Jokinen and A. Mazumdar, JCAP 0610, 007 (2006).
  * (32) R. Allahverdi and A. Mazumdar, JCAP 0610, 008 (2006). M. Berkooz, D. J. H. Chung and T. Volansky, Phys. Rev. Lett. 96, 031303 (2006). R. Allahverdi and A. Mazumdar, arXiv:0802.4430 [hep-ph].
  * (33) A. A. Starobinsky, JETP Lett. 30, 682 (1979) [Pisma Zh. Eksp. Teor. Fiz. 30, 719 (1979)]. S. Dodelson, W. H. Kinney and E. W. Kolb, Phys. Rev. D 56, 3207 (1997); M. Maggiore, Phys. Rept. 331, 283 (2000); T. L. Smith, M. Kamionkowski and A. Cooray, Phys. Rev. D 73, 023504 (2006); arXiv:0802.1530 [astro-ph].
  * (34) M. Kamionkowski, A. Kosowsky and M. S. Turner, Phys. Rev. D 49, 2837 (1994); C. Grojean and G. Servant, Phys. Rev. D 75, 043507 (2007).
  * (35) B. J. Owen and B. S. Sathyaprakash, Phys. Rev. D 60, 022002 (1999).
  * (36) P. de Bernardis, M. Bucher, C. Burigana, L. Piccirillo and f. P. Collaboration, arXiv:0808.1881 [astro-ph].


# Using Joint Random Partition Models for Flexible Change Point Analysis in
Multivariate Processes

José J. Quinlan
Department of Statistics, Pontificia Universidad Católica de Chile, Chile
<EMAIL_ADDRESS>Garritt L. Page
Department of Statistics, Brigham Young University, USA
<EMAIL_ADDRESS>Luis M. Castro
Department of Statistics, Pontificia Universidad Católica de Chile, Chile
Millennium Nucleus Center for the Discovery of Structures in Complex Data,
Chile
Centro de Riesgos y Seguros UC, Pontificia Universidad Católica de, Chile,
Chile
<EMAIL_ADDRESS>
###### Abstract

Change point analyses are concerned with identifying positions of an ordered
stochastic process that undergo abrupt local changes of some underlying
distribution. When multiple processes are observed, it is often the case that
information regarding the change point positions is shared across the
different processes. This work describes a method that takes advantage of this
type of information. Since the number and position of change points can be
described through a partition with contiguous clusters, our approach develops
a joint model for these types of partitions. We describe computational
strategies associated with our approach and illustrate improved performance in
detecting change points through a small simulation study. We then apply our
method to a financial data set of emerging markets in Latin America and
highlight interesting insights discovered due to the correlation between
change point locations among these economies.

Key words: correlated random partitions, multiple change point analysis,
multivariate time series.

##  1 Introduction

Change point analyses identify times or positions of an ordered stochastic
process that undergo abrupt local changes. These abrupt changes are typically
seen as shifts in expectation, variability, or shape of an underlying
distribution (or some combination of the three). Methods that detect change
points have been employed in a variety of fields, including finance (Wood et
al., 2021), climatology (Gupta et al., 2021), and ecology (Jones et al., 2021)
to name a few. Due to this, many change point methods have been proposed in
the statistical literature both in a univariate (see for example Arellano-
Valle et al., 2013) and a multivariate setting (see Truong et al., 2020, for a
comprehensive review).

The phenomenon that motivates our research is the so-called financial
contagion or simply contagion. This phenomenon can be understood as the spread
of financial crises from one country to another (see for example Lowell et
al., 1998; Valdés, 2000; de P. Filleti et al., 2008, among others). To
illustrate contagion consider the price and returns of the five markets
displayed in Figure 1. Note that the overall trend of the price in the Latin
American markets (Argentina, Brazil, Chile, and Mexico) seem to coincide.
However, the USA market (Dow Jones index) presents a different trend during
the same observation period (1995 to 2001). It is important to note that, in
the second half of 1998, Dow Jones suffered a slight crash (a change in
volatility according to Figure 1, left column) due to the Russian financial
crisis and the Long Term Capital Management episode. Consequently, we
hypothesize that a change point in a mature market such as the US could
produce change points in emerging markets such as those from Latin American,
or simply, the financial contagion between the US market and Latin American
markets could increase the chance of change points occurring in the later
markets when they occur in the former. Consequently, the method we develop
will incorporate dependence between change point probabilities across multiple
processes which could potentially improve the ability of detecting a change
point compared to an independent model.

Figure 1: Daily returns and price of stock market indexes of Argentina,
Brazil, Chile, Mexico and USA. Returns are calculated using
$R_{t}=(P_{t}-P_{t-1})/P_{t-1}$ for $t\in\\{2,\ldots,n\\}$, where $P_{t}$ is
the price of day $t$.

One commonly used approach in the statistical literature for detecting time-
series change points is based on product partition models (PPM). These models,
which were introduced by Barry and Hartigan (1992), assume that (a) the number
and positions of change points are random and, (b) observations within the
same block are assumed to follow the same distribution. Thus, the inferential
problem reduces to identifying a partition where each cluster is a collection
of consecutive data points and then estimate parameters associated with each
cluster’s assumed data model. From a Bayesian viewpoint, a prior distribution
on the space of partitions, which are restricted to be contiguous, is needed.
Barry and Hartigan (1992) use a prior for which prior change point
probabilities are based on the so-called cohesion function studied in Yao
(1984). Since Barry and Hartigan (1992) many other PPM type approaches to
change point analysis have be developed (see, for example Loschi and Cruz,
2002; Loschi et al., 2003; Loschi and Cruz, 2005; Loschi et al., 2005, 2010;
Martínez and Mena, 2014; García and Gutiérrez-Peña, 2019; Pedroso et al.,
2021).

Most of the existing approaches based on PPMs for detecting change points in
time series treat the series individually. However, in the presence of
contagion, the information available from several series could improve the
accuracy of the change point detection mechanism compared to when series are
treated individually. In general, the strategies for detecting change points
in the multivariate context focus on detecting changes in the joint
distribution of the coordinates of a multivariate process across time. For
example, Cheon and Kim (2010) developed a Bayesian model for detecting changes
in the mean and variance when data follow a multivariate normal distribution.
Their approach considers a latent vector that identifies change point
positions for partitioning the observations. Nyamundanda et al. (2015)
proposed a Bayesian latent variable PPM, the so-called product partition
latent variable model (PPLVM), providing a flexible framework for detecting
multiple change point detection in multivariate data. The key feature of the
PPLVM is that it can be used to detect distributional changes in the mean and
covariance of the series, even in high-dimensional settings. Recently, Jin et
al. (2021) developed a hierarchical Bayesian model for changes in the mean
process were the prior distribution for change point configurations arises
from a Poisson-Dirichlet process, restricted to contiguous partitions (see
Martínez and Mena, 2014, for more details).

In contrast to the approaches described above, there are models that do not
introduce explicitly a distribution for partitions of contiguous clusters. For
example, Killick et al. (2012) considered a procedure for detecting change
points by minimizing a cost function, which detects the optimal number and
location of change points with a linear computational complexity under mild
conditions. Tveten et al. (2021) also minimize a cost function when searching
for change points in cross-correlated processes. Matteson and James (2014)
proposed a robust nonparametric method using a divergence measure based on
Euclidean distances. With this method, the authors showed that it is possible
to detect any distributional change within an independent sequence of random
variables without making any distributional assumption beyond the existence of
the $\alpha$-th absolute moment, for some $\alpha\in(0,2)$. Padilla et al.
(2021) proposed a novel change point detection algorithm based on the
Kolmogorov–Smirnov statistic and showed that it is nearly minimax rate optimal
under suitable conditions.

Other approaches for detecting changes in a multivariate process that is more
inline with what we proposed assume that each process has its own change point
structure. Harlé et al. (2016) introduce a set of independent binary vectors
whose entries indicate which coordinate of the multivariate process changes.
They do this using a composite marginal likelihood based on Wilcoxon’s rank-
sum test and a suitable prior for the binary vectors. Their approach allows us
to include dependence among change points between different coordinates.
Moreover, Fan and Mackey (2017) introduce a set of change point indicator
variables for each coordinate and time such that the prior change point
probabilities for all coordinates at a fixed time are the same but change
through time.

Our approach, which is motivated by the contagion phenomena, considers
simultaneous changes in all parameters associated with a particular process
dealing with multiple univariate time series (not all of which are necessarily
the same type of response). Thus, each series has its change point structure
where the change points between them may or may not coincide in time. We base
our approach on elements from the method described in Fan and Mackey (2017)
combined with a novel multivariate extension of the PPM approach of Barry and
Hartigan (1992). The resulting method takes advantage of the existing
correlation in change point locations between series. This strategy requires
specifying a joint prior distribution for a collection of partitions.
Constructing these types of dependent partition models over a series of
partitions has only very recently been considered in the literature by, for
instance, Zanini et al. (2019); Page et al. (2021). The method we present is
the first work that we are aware of that considers jointly modeling contiguous
partitions.

The remainder of the article is organized as follows. Section 2 provides
notation and background to the change point PPM. In Section 3 we describe our
approach of incorporating dependence between change point probabilities,
provide some theoretical properties and details with regards to computation.
Section 4 describes a numerical experiment designed to study our method’s
ability to detect change points, and in Section 5 we apply our approach to the
finance data concerning stock market returns of five countries. We close the
paper with some concluding remarks in Section 6.

##  2 Background and Preliminaries

To make the paper self-contained, we start this section with a background
related to PPMs, introducing some notation we will use throughout the paper.

###  2.1 Partition Definition and Notation

Without loss of generality, consider $L>1$ time series
$\bm{y}_{i}=(y_{i,1},\ldots,y_{i,n})^{\top}$, each of length $n>2$. Change
points occur when the behavior of $\bm{y}_{i}$ undergoes sudden changes at
unknown times. These times of sudden changes partition $\\{1,\ldots,n\\}$ into
$k_{i}$ contiguous sets, say $\rho_{i}=\\{S_{i,1},\ldots,S_{i,k_{i}}\\}$, for
some $k_{i}\in\\{1,\ldots,n\\}$. Here, $S_{i,j}$ is the $j$th block and
$k_{i}$ is the number of blocks in $\rho_{i}$. The space of these types of
partitions will be denoted by $\mathcal{C}_{n}$.

An alternative, and commonly used, way to denote a partition $\rho_{i}$ is
through a collection of $n$ cluster labels
$\bm{e}_{i}=(e_{i,1},\ldots,e_{i,n})^{\top}$, such that $e_{i,t}=j$ if $t\in
S_{i,j}$. Due to the contiguous structure of the $\rho_{i}$ that we consider,
the cluster labels that define $\rho_{i}$ necessarily are such that
$1=e_{i,1}\leq e_{i,2}\leq\cdots\leq e_{i,n}=k_{i}$ and
$e_{i,t}-e_{i,t-1}\in\\{0,1\\}$, for all $t\in\\{2,\ldots,n\\}$. Therefore,
change point locations can be identified using $\bm{e}_{i}$. Letting
$\tau_{i,0}=0$ and $\tau_{i,j}=\max\\{t\in\\{1,\ldots,n\\}:e_{i,t}=j\\}$ for
all $j\in\\{1,\ldots,k_{i}\\}$, we have that
$S_{i,j}=\\{\tau_{i,j-1}+1,\ldots,\tau_{i,j}\\}$ and the set
$\bm{\tau}_{i}=\\{\tau_{i,1}+1,\ldots,\tau_{i,k_{i}-1}+1\\}$ identifies the
locations at which change points in $\bm{y}_{i}$ occur.

In a change point setting there is another way to represent $\rho_{i}$ that
will prove to be very useful for our approach. It is based on a set of change
point indicators $\bm{c}_{i}=(c_{i,1},\ldots,c_{i,n-1})^{\top}$, such that
$c_{i,t}=1$ if time $t+1$ is a change point in $\bm{y}_{i}$, and $c_{i,t}=0$
otherwise. Note that $\bm{e}_{i}$ and $\bm{c}_{i}$ have a one-to-one
correspondence through the relation $c_{i,t}=e_{i,t+1}-e_{i,t}$, for all
$t\in\\{1,\ldots,n-1\\}$, where $e_{i,1}=1$ by construction. The number of
change points can be identified using $\bm{c}_{i}$ by noticing that
$k_{i}=1+\sum_{t=1}^{n-1}c_{i,t}$. In what follows, we will use $\rho_{i}$,
$\bm{e}_{i}$ and $\bm{c}_{i}$ interchangeably. With the necessary notation
introduced, we next describe the change point PPM.

###  2.2 Change Point Product Partition Models

For the $i$-th sequence, the PPM is a discrete distribution on space
$\mathcal{C}_{n}$ such that

$\mathbb{P}(\rho_{i}=\\{S_{i,1},\ldots,S_{i,k_{i}}\\})=\frac{\prod_{j=1}^{k_{i}}c(S_{i,j})}{\sum_{\rho_{i}\in\mathcal{C}_{n}}\prod_{j=1}^{k_{i}}c(S_{i,j})},$

where $c(S_{i,j})$ is referred to as a cohesion function and measures the a
priori belief that elements in $S_{i,j}$ co-cluster. The change point PPM as
described in Barry and Hartigan (1992), Loschi et al. (2003), and others, uses
Yao (1984)’s cohesion function to assign probabilities to each element in
$\mathcal{C}_{n}$. Yao (1984)’s cohesion function applied to contiguous
$S_{i,j}$ results in

$c(S_{i,j}\,;\,p_{i})=\left\\{\begin{array}[]{lc}p_{i}(1-p_{i})^{\tau_{i,j}-\tau_{i,j-1}-1},&\mathrm{if}\quad\tau_{i,j}<n\\\
(1-p_{i})^{\tau_{i,j}-\tau_{i,j-1}-1},&\mathrm{if}\quad\tau_{i,j}=n\end{array}\right.,$
(1)

for some $p_{i}\in[0,1]$ such that $p_{i}=\mathbb{P}(c_{i,t}=1\mid p_{i})$.
Based on this cohesion we have that

$\sum_{\rho_{i}\in\mathcal{C}_{n}}\prod_{j=1}^{k_{i}}c(S_{i,j}\,;\,p_{i})=1.$

Thus, the change point PPM takes on the following form

$\displaystyle\mathbb{P}(\rho_{i}=\\{S_{i,1},\ldots,S_{i,k_{i}}\\}\mid
p_{i})=\prod_{j=1}^{k_{i}}c(S_{i,j}\,;\,p_{i})=p_{i}^{k_{i}-1}(1-p_{i})^{n-k_{i}}.$

Once the partition model has been specified, the key idea behind change point
modeling from a partition perspective is that observations within the same
block are assumed to follow a common distribution, whereas different
distributions are assumed between blocks. Following Barry and Hartigan
(1992)’s approach, given $\rho_{i}$, the joint density of $\bm{y}_{i}$ is
written as a product of $k_{i}$ data factors, also known as marginal
likelihoods, which measure the similarity of observations within each block.
More precisely,

$\displaystyle\begin{split}f(\bm{y}_{i}\mid\rho_{i},\bm{\xi}_{i})&=\prod_{j=1}^{k_{i}}\mathcal{F}_{i}(\bm{y}_{i,j}\mid\bm{\xi}_{i}),\\\
\mathcal{F}_{i}(\bm{y}_{i,j}\mid\bm{\xi}_{i})&=\int_{\bm{\Theta}_{i}}\mathcal{L}_{i}(\bm{y}_{i,j}\mid\bm{\theta}_{i},\bm{\xi}_{i})\mathrm{d}G_{i}(\bm{\theta}_{i}\mid\bm{\xi}_{i}),\end{split}$
(2)

where $\bm{y}_{i,j}=(y_{i,t}:t\in S_{i,j})^{\top}$ and
$\mathcal{L}_{i}(\,\,\cdot\mid\bm{\theta}_{i},\bm{\xi}_{i})$ is a likelihood
function indexed by the set of parameters $\bm{\theta}_{i}\in\bm{\Theta}_{i}$
which are block-specific, and $\bm{\xi}_{i}\in\bm{\Xi}_{i}$ a collection of
parameters that are common to all blocks. Further,
$G_{i}(\,\,\cdot\mid\bm{\xi}_{i})$ is a suitable prior for $\bm{\theta}_{i}$.
The data generating mechanism (2) along with prior distributions for
$\rho_{i}$ and $\bm{\xi}_{i}$ (if applicable) completely specify the Bayesian
change point PPM.

It is common to select
$\mathcal{L}_{i}(\,\,\cdot\mid\bm{\theta}_{i},\bm{\xi}_{i})$ and
$G_{i}(\,\,\cdot\mid\bm{\xi}_{i})$ such that they form a conjugate pair which
results in $\mathcal{F}_{i}(\bm{y}_{i,j}\mid\bm{\xi}_{i})$ being available in
closed form. That said, the choices for
$\mathcal{L}_{i}(\,\,\cdot\mid\bm{\theta}_{i},\bm{\xi}_{i})$ and
$G_{i}(\,\,\cdot\mid\bm{\xi}_{i})$ in (2) can be quite general, depending on
the nature of $\bm{y}_{i}$ and $\bm{\theta}_{i}$. Examples of this are data
that follow an Ornstein-Uhlenbeck process with a Normal-Gamma prior for mean-
precision parameters (Martínez and Mena, 2014) (here, $\bm{\xi}_{i}$ is a case
dependency parameter with a $\mathrm{Uniform}(0,1)$ prior) and independent
data belonging to the exponential family with a conjugate prior for the
natural parameters (Loschi and Cruz, 2005) (in this case, there is no
$\bm{\xi}_{i}$). The types of marginal likelihoods just described (and others
not listed) are easily applied using our method. Even so, in what follows, we
will focus on the following specification (which is suitable for changes in
mean and variance for data supported on $\mathbb{R}$). Let
$\bm{\theta}_{i}=(\mu_{i},\sigma_{i}^{2})^{\top}\in\bm{\Theta}_{i}=\mathbb{R}\times(0,+\infty)$
and

$\displaystyle\begin{split}\mathcal{L}_{i}(\bm{y}_{i,j}\mid\bm{\theta}_{i})&=\prod_{t\in
S_{i,j}}\mathrm{Normal}(y_{i,t}\mid\mu_{i},\sigma_{i}^{2}),\\\
G_{i}(\bm{\theta}_{i})&=\mathrm{Normal}(\mu_{i}\mid\mu_{i,0},(\kappa_{i,0})^{-1}\sigma_{i}^{2})\mathrm{Inv{\textrm{-}}Gamma}(\sigma_{i}^{2}\mid\alpha_{i,0},\beta_{i,0}).\end{split}$
(3)

Here, $\mu_{i,0}\in\mathbb{R}$ and $\kappa_{i,0},\alpha_{i,0},\beta_{i,0}>0$
are fixed hyperparameters, $\mathrm{Normal}(\,\,\cdot\mid\mu,\sigma^{2})$
denotes a normal density with mean $\mu\in\mathbb{R}$ and variance
$\sigma^{2}>0$, and $\mathrm{Inv{\textrm{-}}Gamma}(\,\,\cdot\mid\alpha,\beta)$
denotes an inverse Gamma density with shape $\alpha>0$ and scale $\beta>0$. It
is well known that the $j$th data factor induced by (3) is given by

$\mathcal{F}_{i}(\bm{y}_{i,j})=t_{n_{i,j}}(\bm{y}_{i,j}\mid
2\alpha_{i,0},\mu_{i,0}\mathbf{1}_{n_{i,j}},(\alpha_{i,0})^{-1}\beta_{i,0}\\{\mathbf{I}_{n_{i,j}}+(\kappa_{i,0})^{-1}\mathbf{J}_{n_{i,j}}\\}),$
(4)

where $n_{i,j}$ is the cardinality of $S_{i,j}$, $\bm{1}_{p}\in\mathbb{R}^{p}$
is the vector with entries equal 1, $\mathbf{I}_{p}\in\mathbb{R}^{p\times p}$
is the identity matrix and $\mathbf{J}_{p}\in\mathbb{R}^{p\times p}$ is the
matrix with all entries equal 1. Also,
$t_{p}(\,\,\cdot\mid\nu,\bm{\mu},\bm{\Sigma})$ is the $p$-dimensional
Student’s $t$-density with degrees of freedom $\nu>0$, location vector
$\bm{\mu}\in\mathbb{R}^{p}$ and scale matrix
$\bm{\Sigma}\in\mathbb{S}^{p\times p}$, where $\mathbb{S}^{p\times p}$ denotes
the space of positive-definite matrices.

The set of hyperparameters $(\mu_{i,0}$, $\kappa_{i,0}$,$\alpha_{i,0}$,
$\beta_{i,0})^{\top}$ play a crucial role in determining what constitutes a
change point. For example, setting $\alpha_{i,0}$ close to one will result in
(4) having thick tails so that a change point would necessarily need to be far
from the center. Conversely, with a large value of $\alpha_{i,0}$ (4)
approximates a normal distribution and points not far from the center can
still be change points. Consequently, thought must be dedicated to assigning
values to the marginal likelihood parameters. In Section 3.2 we discuss an
empirical Bayes method that produces reasonable values for them in absence of
prior information.

##  3 The Joint Prior Distribution on a Collection of Partitions

We now describe our approach of formulating a joint model for a sequence of
partitions. As mentioned, the partition $\rho_{i}$ for the $i$-th time series
$\bm{y}_{i}$ has a one-to-one correspondence with $\bm{c}_{i}$. Thus, any
prior distribution for $\mathbf{C}=(\bm{c}_{1},\ldots,\bm{c}_{L})$, say
$\pi(\mathbf{C})$, uniquely determines a prior for
$\bm{\rho}=(\rho_{1},\ldots,\rho_{L})$. We start by describing our joint model
as an extension of the change point PPM and then we connect it to
$\pi(\mathbf{C})$, which is what we ultimately use in our approach as it
facilitates computation.

In our setup, rather than consider a single probability parameter $p_{i}$ for
the $i$-th series, we define
$\bm{\tilde{p}}_{i}=(p_{i,1},\ldots,p_{i,n-1})^{\top}$ with $p_{i,t}\in[0,1]$
and extend the cohesion in (1) to

$c^{\star}(S_{i,j}\,;\,\bm{\tilde{p}}_{i})=\left\\{\begin{array}[]{lc}p_{i,\tau_{i,j}}\prod_{t=\tau_{i,j-1}+1}^{\tau_{i,j}-1}(1-p_{i,t}),&\mathrm{if}\quad\tau_{i,j}<n,\\\
\prod_{t=\tau_{i,j-1}+1}^{\tau_{i,j}-1}(1-p_{i,t}),&\mathrm{if}\quad\tau_{i,j}=n.\end{array}\right.$
(5)

Using the cohesion (5) for contiguous partitions still results in
$\sum_{\rho_{i}\in\mathcal{C}_{n}}\prod_{j=1}^{k_{i}}c^{\star}(S_{i,j}\,;\,\bm{\tilde{p}}_{i})=1.$
Therefore, the partition probabilities become

$\mathbb{P}(\rho_{i}=\\{S_{i,1},\ldots,S_{i,k_{i}}\\}\mid\bm{\tilde{p}}_{i})=\prod_{j=1}^{k_{i}}c^{\star}(S_{i,j}\,;\,\bm{\tilde{p}}_{i})=\prod_{t\in\bm{T}_{i}}p_{i,t}\prod_{t\notin\bm{T}_{i}}(1-p_{i,t}),$

where $\bm{T}_{i}=\\{\tau_{i,1},\ldots,\tau_{i,k_{i}-1}\\}$. Including all $L$
partitions, the joint partition model becomes

$\displaystyle\mathbb{P}(\rho_{i}=\\{S_{i,1},\ldots,S_{i,k_{i}}\\}:i=1,\ldots,L\mid\bm{\tilde{p}}_{1},\ldots,\bm{\tilde{p}}_{L})$
$\displaystyle=\prod_{i=1}^{L}\mathbb{P}(\rho_{i}=\\{S_{i,1},\ldots,S_{i,k_{i}}\\}\mid\bm{\tilde{p}}_{i})$
$\displaystyle=\prod_{i=1}^{L}\left\\{\prod_{t\in\bm{T}_{i}}p_{i,t}\prod_{t\notin\bm{T}_{i}}(1-p_{i,t})\right\\}$

It is important to stress that our extension (i.e., introducing $p_{i,t}$)
provides more flexibility for modeling simultaneous change point
configurations $\bm{\rho}$ by allowing us to correlate probabilities of a
change point at each time point. As a consequence, the change point indicators
$c_{i,t}$ are assumed conditionally independent with their own probability
$p_{i,t}\in(0,1)$ of detecting a change $(c_{i,t}=1)$.

Let $\bm{p}_{t}=(p_{1,t},\ldots,p_{L,t})^{\top}$ be a $L$-dimensional vector
of probabilities supported on the space $(0,1)^{L}$. To specify a multivariate
distribution for $\bm{p}_{t}$, we consider the bijective transformation
$\mathrm{logit}(\bm{p}_{t})=\big{(}\log\big{(}\frac{p_{1,t}}{1-p_{1,t}}\big{)},\ldots,\log\big{(}\frac{p_{L,t}}{1-p_{L,t}}\big{)}\big{)}^{\top}$
which is defined on the Euclidean space $\mathbb{R}^{L}$, and model it with a
multivariate Student’s-$t$ distribution. The reason for selecting a
Student’s-$t$ distribution instead of, for example, a multivariate normal
distribution is that extreme probabilities (near 0 or 1) are more achievable
due to the thicker tails of the Student’s-$t$. In summary, the proposed model
for partitions $\bm{\rho}$ can be formulated using the following hierarchical
structure:

$\displaystyle\begin{split}c_{i,t}\mid
p_{i,t}&\stackrel{{\scriptstyle\mathrm{ind.}}}{{\sim}}p_{i,t}^{c_{i,t}}(1-p_{i,t})^{1-c_{i,t}},\\\
\mathrm{logit}(\bm{p}_{t})&\stackrel{{\scriptstyle\mathrm{iid.}}}{{\sim}}t_{L}(\nu_{0},\bm{\mu}_{0},\bm{\Sigma}_{0}),\end{split}$
(6)

with $\bm{\mu}_{0}=(\mu^{0}_{1},\ldots,\mu^{0}_{L})^{\top}$,
$\bm{\Sigma}_{0}=(\sigma^{0}_{l_{1},l_{2}}:l_{1},l_{2}\in\\{1,\ldots,L\\}).$
In what follows, we will refer to the model comprised of (3), (5), and (6) as
the correlated change point product partition model or simply, CCP-PPM.

Note that, specifying adequate values for $\nu_{0}$, $\bm{\mu}_{0}$ and
$\bm{\Sigma}_{0}$ in (6) must be done with caution as $p_{i,t}$ are not
invariant to their selection. In the absence of information regarding these
parameters, we provide an interesting approach for selecting them in Section
3.2.

###  3.1 Properties of the Joint Model on Partitions

In this section, we provide some interesting properties that are consequences
of modeling the change point indicators with (6). The proofs of all
propositions are provided in the Appendix.

######  Proposition 1.

Under the assumptions in (6), each of the $\bm{c}_{i},$ $i=1,\ldots,L,$
follows a change point PPM based on Yao’s cohesion with probability parameter

$\phi_{i}=\int_{\mathbb{R}}\left\\{\frac{\exp(z)}{1+\exp(z)}\right\\}t_{1}(z\mid\nu_{0},\mu^{0}_{i},\sigma^{0}_{i,i})\mathrm{d}z.$
(7)

A consequence of Proposition 1 is that the correlation in the change point
probabilities from our model only exists across the $L$ series for a fixed
$t$. Within a series, the probability of a change point at time $t_{1}$ is
independent of time $t_{2}$.

The next proposition provides an interesting result regarding the number of
expected change points based on the CCP-PPM.

######  Proposition 2.

Under the assumptions in (6), the number of change points $(k_{i}-1)$ for the
$i$th series satisfies $(k_{i}-1)\sim\mathrm{Binomial}(n-1,\phi_{i})$, where
$\phi_{i}$ is given by (7). Additionally, for series $i$ and $s$, the
distribution for $(k_{i}-1,k_{s}-1)^{\top}$ is a mixture of a product of two
Poisson’s Binomial distributions (Wang, 1993) with
$\mathrm{Cov}(k_{i}-1,k_{s}-1)=(n-1)(\varphi_{i,s}-\phi_{i}\phi_{s})$, where

$\varphi_{i,s}=\int_{\mathbb{R}^{2}}\left[\frac{\exp(z_{i})\exp(z_{s})}{\\{1+\exp(z_{i})\\}\\{1+\exp(z_{s})\\}}\right]t_{2}(\bm{z}_{A}\mid\nu_{0},\bm{\mu}_{A,0},\bm{\Sigma}_{A,0})\mathrm{d}\bm{z}_{A}.$
(8)

Here, $\bm{z}_{A}=(z_{i},z_{s})^{\top}$,
$\bm{\mu}_{A,0}=(\mu^{0}_{i},\mu^{0}_{s})^{\top}$ and
$\bm{\Sigma}_{A,0}=(\sigma^{0}_{l_{1},l_{2}}:l_{1},l_{2}\in\\{i,s\\})$.

Note that $\mathrm{E}(k_{i}-1)=(n-1)\phi_{i}$ easily follows from Proposition
2. For the values of $\mu^{0}_{i}$ and $\sigma^{0}_{i,i}$ considered in Figure
2, the expected number of change points in each series ranges between 6.19 and
12.29 a priori. The bottom row of Figure 2 provides values for
$\mathrm{Corr}(k_{i}-1,k_{s}-1)$ based on a numerical approximation of
$\varphi_{i,s}$. As expected the number of change points in two series with
high correlation between $\mathrm{logit}(p_{i,t})$ and
$\mathrm{logit}(p_{s,t})$ will be similar.

The last proposition derives the conditional probabilities of change point
indicators. These are of particular interest as they illustrate how the
probability of a change point across series varies as one series experiences a
change point.

######  Proposition 3.

Under the assumptions in (6), the probability of a change point occurring at
time $t$ in the $i$th series given that one occurred at time $t$ in the $s$th
series is

$\mathbb{P}(c_{i,t}=1\mid c_{s,t}=1)=\frac{\varphi_{i,s}}{\phi_{s}}.$ (9)

Here $\phi_{s}$ and $\varphi_{i,s}$ are given by (7) and (8), respectively.

The top row of Figure 2 provides values for $\mathbb{P}(c_{i,t}=1\mid
c_{s,t}=1)$. The integral in $\mathbb{P}(c_{i,t}=1\mid c_{s,t}=1)$ was
approximated using the statistical software R (R Core Team 2021). As expected,
the higher the correlation between $\mathrm{logit}(p_{i,t})$ and
$\mathrm{logit}(p_{s,t})$ the higher the conditional probabilities a priori,
for the values of $\mu^{0}_{i}$ and $\sigma^{0}_{i,i}$ considered.

Figure 2: The top row displays values of 9 for a few values of $\mu_{i}^{0}$
and $\sigma_{i}^{0}$ and in increasing sequence of correlations between
$\mathrm{logit}(p_{i,t})$ and $\mathrm{logit}(p_{s,t})$. The bottom row
displays $\mathrm{Corr}(k_{i}-1,k_{s}-1)$ for the same values of
$\mu_{i}^{0}$, $\sigma_{i}^{0}$, and correlations between
$\mathrm{logit}(p_{i,t})$ and $\mathrm{logit}(p_{s,t})$. The integral in (8)
was approximated numerically.

###  3.2 Selection of tuning parameters

Like most change point methods, the posterior probability of classifying a
point as a change point can be sensitive to marginal likelihood specification
and prior parameter selection. We refer to these parameters as tuning
parameters. In some cases, the practitioner can inform the procedure regarding
a change point, which guides tuning parameter selection. Without this
information, it is appealing to have a procedure that produces default values
for the tuning parameters. Therefore, we describe an empirical Bayes approach
to selecting values for $(\mu_{i,0}$, $\kappa_{i,0}$, $\alpha_{i,0}$,
$\beta_{i,0})^{\top}$. The approach we describe is geared towards situations
in which the magnitudes of change points are relatively small.

First note that from (4) we have

$0\leq\mathrm{Corr}(y_{i,r},y_{i,s})=\frac{\mathbb{P}(e_{i,r}=e_{i,s})}{1+\kappa_{i,0}}\leq\frac{1}{1+\kappa_{i,0}},$

for all $r,s\in\\{1,\ldots,n\\}$ and $r\neq s$. Under the scenario of no
change points, $\mathrm{Corr}(y_{i,r},y_{i,s})$ is equal to the upper bound
$(1+\kappa_{i,0})^{-1}$, which is constant as a function of $(r,s)$. Thus, a
value for $\kappa_{i,0}$ could be empirically selected using
$\mathrm{Corr}(y_{i,r},y_{i,s})$. This correlation can be estimated by
inspecting the empirical $\ell$-lag autocorrelation for $\bm{y}_{i}$, say
$\widehat{c_{i,\ell}}$, and choosing the smallest $\ell\geq 1$ such that
$\widehat{c_{i,\ell}}>0$. Then set
$\kappa_{i,0}=(\widehat{c_{i,\ell}})^{-1}(1-\widehat{c_{i,\ell}})$. From
there, the maximum likelihood estimators of the parameters in the Student’s
$t$-distribution given in (4) based on the entire $\bm{y}_{i}$ can be used to
provide values for $\mu_{i,0}$, $\alpha_{i,0}$, and $\beta_{i,0}$. Let
$d=2\alpha_{i,0}$, $m=\mu_{i,0}$ and
$s=(\alpha_{i,0}\kappa_{i,0})^{-1/2}\\{\beta_{i,0}(\kappa_{i,0}+1)\\}^{1/2}$.
Then, we set $\mu_{i,0}=\widehat{m_{i}}$, $\alpha_{i,0}=0.5\widehat{d_{i}}$
and
$\beta_{i,0}=0.5\widehat{d_{i}}(1-\widehat{c_{i,\ell}})\widehat{s_{i}}^{2}$
where $(\widehat{d_{i}},\widehat{m_{i}},\widehat{s_{i}})^{\top}$ denote the
maximum likelihood estimates of $(d,m,s)^{\top}$.

Now, we focus on $(\nu_{0},\bm{\mu}_{0},\bm{\Sigma}_{0})^{\top}$. Although
$\mathrm{E}(\bm{p}_{t})$ and $\mathrm{Var}(\bm{p}_{t})$ do not exist in a
closed form, a first-order Taylor expansion provides approximations to them.
If $\nu_{0}>2$, then

$\displaystyle\mathrm{E}(\bm{p}_{t})$ $\displaystyle\approx$
$\displaystyle\left(\frac{\exp(\mu^{0}_{1})}{1+\exp(\mu^{0}_{1})},\ldots,\frac{\exp(\mu^{0}_{L})}{1+\exp(\mu^{0}_{L})}\right)^{\top},$
$\displaystyle\mathrm{Var}(\bm{p}_{t})$ $\displaystyle\approx$
$\displaystyle\left(\frac{\nu_{0}}{\nu_{0}-2}\right)\mathbf{J}(\bm{\mu}_{0})\bm{\Sigma}_{0}\mathbf{J}(\bm{\mu}_{0}),$
$\displaystyle\mathbf{J}(\bm{\mu}_{0})$ $\displaystyle=$
$\displaystyle\mathrm{diag}\left(\frac{\exp(\mu^{0}_{1})}{\\{1+\exp(\mu^{0}_{1})\\}^{2}},\ldots,\frac{\exp(\mu^{0}_{L})}{\\{1+\exp(\mu^{0}_{L})\\}^{2}}\right).$

After choosing prior guesses for $\mathrm{E}(\bm{p}_{t})$ and
$\mathrm{Var}(\bm{p}_{t})$, say
$\bm{m}_{0}=(m^{0}_{1},\ldots,m^{0}_{L})^{\top}$ and $\mathbf{S}_{0}$
respectively, we set

$\displaystyle\bm{\mu}_{0}$ $\displaystyle=$
$\displaystyle\left(\log\left(\frac{m^{0}_{1}}{1-m^{0}_{1}}\right),\ldots,\log\left(\frac{m^{0}_{L}}{1-m^{0}_{L}}\right)\right)^{\top},$
$\displaystyle\bm{\Sigma}_{0}$ $\displaystyle=$
$\displaystyle\left(\frac{\nu_{0}-2}{\nu_{0}}\right)\mathbf{D}(\bm{m}_{0})^{-1}\mathbf{S}_{0}\mathbf{D}(\bm{m}_{0})^{-1},$
$\displaystyle\mathbf{D}(\bm{m}_{0})^{-1}$ $\displaystyle=$
$\displaystyle\mathrm{diag}\left(\frac{1}{m^{0}_{1}(1-m^{0}_{1})},\ldots,\frac{1}{m^{0}_{L}(1-m^{0}_{L})}\right).$

We recommend setting $\nu_{0}=3$, which is the least integer such that the
approximations described above exist. In the case that no prior information is
available to guide specifying $\bm{m}_{0}$ and $\mathbf{S}_{0}$, the following
empirical approach can be used. Set $m^{0}_{1}=\cdots=m^{0}_{L}=n^{-1}$. For
$\mathbf{S}_{0}$, a compound symmetry covariance matrix
$\sigma^{2}_{0}\\{(1-r_{0})\mathbf{I}_{L}+r_{0}\mathbf{J}_{L}\\}$ can be used,
where $\sigma^{2}_{0}=m^{0}_{1}(1-m^{0}_{1})/n=n^{-3}(n-1)$ and $r_{0}=0.5$.

###  3.3 Posterior Sampling

The joint posterior distribution of $\bm{p}_{t}$ and $\mathbf{C}$ is not
analytically tractable. Therefore we resort to sampling from it using Markov
Chain Monte Carlo (MCMC) methods. The MCMC algorithm we construct is very
straightforward to implement and is an extension of the Gibbs sampler
described in Loschi et al. (2003) with the main difference being that we
employ Metropolis updates when conditional conjugacy is not available. Before
model fitting, we recommend scaling each series to have mean zero and standard
deviation one. The full conditional distributions we need in our algorithm are
described next.

  * •

For each $i\in\\{1,\ldots,L\\}$ and $t\in\\{1,\ldots,n-1\\}$ update $p_{i,t}$
according to its full conditional density $\pi(p_{i,t}\mid\cdots)$, which is
proportional to

$\frac{1}{p_{i,t}}\left(\frac{p_{i,t}}{1-p_{i,t}}\right)^{c_{i,t}}\left\\{1+\frac{(\mathrm{logit}(\bm{p}_{t})-\bm{\mu}_{0})^{\prime}\bm{\Sigma}_{0}^{-1}(\mathrm{logit}(\bm{p}_{t})-\bm{\mu}_{0})}{\nu_{0}}\right\\}^{-\frac{\nu_{0}+L}{2}}1(p_{i,t}\in(0,1)).$

Here, $1(\,\,\cdot\,\in S)$ is the indicator function of the set $S$. To
update $p_{i,t}$, we employ a random walk Metropolis step with a normal
centered at the previous iteration’s value as a candidate density. The
standard deviation of the normal candidate density is set to 0.005 which
produces an acceptance rate in the general range of 0.2 and 0.5

  * •

For each $i\in\\{1,\ldots,L\\}$, $t\in\\{1,\ldots,n-1\\}$ and $a\in\\{0,1\\}$,
define the set of change point indicators
$\bm{c}^{\scriptscriptstyle{(a)}}=(c_{1}^{\scriptscriptstyle{(a)}},\ldots,c_{n-1}^{\scriptscriptstyle{(a)}})^{\top}$
such that

$c_{s}^{\scriptscriptstyle{(a)}}=\left\\{\begin{array}[]{ll}c_{i,s},&\mbox{if
$s\neq t$}\\\ a,&\mbox{if $s=t$.}\end{array}\right.$

Using $\bm{c}^{\scriptscriptstyle{(a)}}$, we construct the corresponding set
of cluster labels
$\bm{e}^{\scriptscriptstyle{(a)}}=(e_{1}^{\scriptscriptstyle{(a)}},\ldots,e_{n}^{\scriptscriptstyle{(a)}})^{\top}$.
Then, after computing

$\varpi_{i,t}=\frac{\mathbb{P}(c_{i,t}=1\mid\cdots)}{\mathbb{P}(c_{i,t}=0\mid\cdots)}=\frac{\mathcal{F}_{i}(e_{t}^{\scriptscriptstyle{(1)}}\mid\bm{\xi}_{i})\mathcal{F}_{i}(e_{t+1}^{\scriptscriptstyle{(1)}}\mid\bm{\xi}_{i})}{\mathcal{F}_{i}(e_{t}^{\scriptscriptstyle{(0)}}\mid\bm{\xi}_{i})}\left(\frac{p_{i,t}}{1-p_{i,t}}\right),$

where
$\mathcal{F}_{i}(j\mid\bm{\xi}_{i})=\mathcal{F}_{i}(\bm{y}_{i,j}\mid\bm{\xi}_{i})$,
$c_{i,t}$ can be updated using a Bernoulli distribution with probability
parameter

$\mathbb{P}(c_{i,t}=1\mid\cdots)=\frac{\varpi_{i,t}}{1+\varpi_{i,t}}.$

Now, an MCMC algorithm can be obtained by cycling through each of the full
conditionals individually. If a model is proposed so that $\bm{\xi}_{i}$ is
available, it is relatively straightforward to update $\bm{\xi}_{i}$ in the
Gibbs sampler using a Metropolis step. The update is based on the following
full conditional of $\bm{\xi}_{i}$ for each $i\in\\{1,\ldots,L\\}$

$\pi(\bm{\xi}_{i}\mid\cdots)\propto\left\\{\prod_{j=1}^{k_{i}}\mathcal{F}_{i}(\bm{y}_{i,j}\mid\bm{\xi}_{i})\right\\}h_{i}(\bm{\xi}_{i}),$

where $h_{i}$ is a prior density for $\bm{\xi}_{i}$.

##  4 Simulation Study

We conduct a small numerical experiment to study the CCP-PPM’s ability to
detect multiple change points. The experiment is based on generating data sets
containing change points whose times are dependent across series, mimicking
the contagion idea. We consider change points that result from simultaneous
changes in the mean and variance of a normal distribution. Data sets are
generated using three scenarios, with each one producing data sets containing
$L=2$ series of $n=100$ observations. The scenarios are detailed next.

The first scenario referred as data type 1 uses the CCP-PPM as a data
generating mechanism. We set $\nu_{0}=3$, $\bm{\mu}_{0}=(-6,-6)^{\top}$, and
$\bm{\Sigma}_{0}$ to a compound symmetric matrix with variance 10 and
correlation 0.9. Based on these values, (6) is used to create partitions. Once
partitions are formed, we use a $\mathrm{Normal}(0,1)$ to generate cluster
specific means for one series and a $\mathrm{Normal}(4,1)$ for the other.
Cluster specific variances were all generated using an
$\mathrm{Inv{\textrm{-}}Gamma}(10,1)$ distribution. Once cluster specific
parameters were generated, observations were generated using a
$\mathrm{Normal}(0,1)$. An example of this type of data is displayed in the
top row of Figure 3.

Figure 3: Example of the types of data sets used in the simulation study. Top
row corresponds to data type 1, while the second and third rows to data type 2
and data type 3 respectively. Points that are labeled as 1 come from the first
series with those labeled 2 come from the second. The different colors
identify the clusters formed by change points.

The following two scenarios set change point locations in both series at 25,
50, and 75. As a result, the change points of the two series are highly
dependent. Under this setting, four clusters of 25 observations for both
series are obtained. Given this type of partition, we produce observations in
two ways. The first one, which we refer to as data type 2, generates
observations using a normal distribution with the following cluster-specific
means and variances:

  1. -

$\bm{\mu}^{\star}=(-1,0,1,2)^{\top}$ and
$\bm{\sigma}^{\star}=(0.1,0.25,0.5,0.75)^{\top}$,

  2. -

$\bm{\mu}^{\star}=(2,1,0,-1)^{\top}$ and
$\bm{\sigma}^{\star}=(0.1,0.25,0.5,0.75)^{\top}$.

The second scenario, which we refer to as data type 3, generates observations
using a normal distribution and cluster-specific means and variances given by:

  1. -

$\bm{\mu}^{\star}=(-0.25,0,0.25,0.5)^{\top}$ and
$\bm{\sigma}^{\star}=(0.1,0.25,1,0.25)^{\top}$,

  2. -

$\bm{\mu}^{\star}=(-0.25,0,0.25,0.5)^{\top}$ and
$\bm{\sigma}^{\star}=(0.1,2,0.5,1)^{\top}$.

This scenario is included because it provides insight to how the CCP-PPM
approach performs for data similar to that which we consider in Section 5.
Examples of synthetic data sets created from these two scenarios are provided
in the second and third rows of Figure 3. It is important to note that in the
scenario data type 2, the means are the primary driver of change points, while
in the scenario data type 3 the variances, or volatility, is the primary
driver of change points.

We simulated one hundred data sets for each scenario. Then, we fit the CCP-PPM
by collecting 2000 MCMC iterates after discarding the initial 10,000 as burn-
in and thinning by 10 all of which was carried out using the ppmSuite (Page
and Quinlan 2022) R package that is available on CRAN. In addition, we also
consider the following controls for detecting change points:

  * •

The PPM-based method of Wang and Emerson (2015). Since this approach does not
analyze multiple series, we fit it to each series independently instead. The
bcp package (Erdman and Emerson, 2007) found in the R statistical software (R
Core Team, 2021) is used to implement this method. We referred to this method
as the Wang method.

  * •

The method of Matteson and James (2014). For this method, the two series are
analyzed jointly. This method is implemented using the R package ecp (James et
al., 2019). We referred to this method as the Matteson method.

  * •

The method developed in Barry and Hartigan (1992) and Loschi and Cruz (2002).
This method is our most natural competitor and is implemented using the
ppmSuite R package. This method is referred to as the Loschi method. As for
the Wang method, the Loschi method is fit to each series independently.

The tuning parameters used for the CCP-PPM are the following. For the marginal
likelihood we use
$(\mu_{i,0},\kappa_{i,0},\alpha_{i,0},\beta_{i,0})^{\top}=(0,1,2,1)^{\top}$.
These same values are used for the Loschi method. For the prior on
$\bm{p}_{t}$, we use the same values that were used to generate the data type
1. For the Loschi method there is a single $p_{i}$ for each series and we used
$p_{i}\stackrel{{\scriptstyle\mathrm{iid.}}}{{\sim}}\mbox{Beta}(1,20)$.
Default values of the bcp and ecp R-packages were considered for the Wang and
Matteson methods, respectively.

Figure 4: Results from the simulation study. Each boxplot displays the results
for each method based on the 100 datasets generated. For the adjusted Rand
index, sensitivity (true positive rate), and specificity (true negative rate)
higher values indicate superior performance.

After considering the methods described above in each scenario, we classify
any point as a change if its posterior probability of being a change is
greater than 0.5. Then, we compute the overall misclassification rate,
sensitivity (true positive rate), and specificity (true negative rate) of
classifying points as change points. These metrics give us information on the
accuracy of identifying points as change points. In addition, we calculate the
adjusted Rand index (Rand, 1971) (ARI) between the estimated partition based
on change points and the true partition. This metric illustrates how well each
method does at recovering the true partition of contiguous clusters.

The results of the simulation study are provided in Figure 4. Note that the
CCP-PPM outperforms all the considered methods at recovering the true
partition according to the adjusted Rand index regardless of the type of data.
In terms of the overall misclassification rate, our method performs better
than the Loschi and Matteson methods. The Wang method performs comparably way
to our approach in the scenario 1 (data type 1), but in scenario 3 (data type
3) this method performs poorly. Regarding the sensitivity, our approach
outperforms the other three methods in the first and second scenarios (data
type 1 and data type 2, respectively). Nonetheless, in the third scenario
(data type 3) Matteson has the highest sensitivity but at the cost of very low
specificity. In summary, it appears that the CCP-PPM overall performs the best
at detecting the correct number and location of change points.

##  5 Finance Data Application

We now turn our attention to the application that motivated our proposal. As
is commonly done in financial applications, we analyze returns rather than
prices. Returns are defined as $R_{t}=(P_{t}-P_{t-1})/P_{t-1}$ for
$t\in\\{2,\ldots,n\\}$, where $P_{t}$ is the daily price. When analyzing the
data set, we consider contagion both between mature markets and emerging ones
and also between emerging markets. We analyzed the most important Latin
American markets (emerging markets), namely Argentinean, Brazilian, Chilean
and Mexican markets, including the USA market in the analysis (mature market).
Consequently, we fit the CCP-PPM and the change point PPM of Loschi or simply
Loschi method (which is perhaps our method’s most natural competitor, treating
each series independently).

We considered the return series of their main stock indexes, namely, the
MERVAL (Índice de Mercado de Valores de Buenos Aires) of Argentina, the
IBOVESPA (Índice da Bolsa de Valores do Estado de São Paulo) of Brazil, the
IPSA (Índice de Precios Selectivos de Acciones) of Chile, the IPyC (Índice de
Precios y Cotizaciones) of Mexico, and the Dow Jones (Dow Jones Industrial
Average) of USA. The stock returns were recorded daily from October 31, 1995
to October 31, 2000.

We employ the procedure described in Section 3.2 to produce values for the
tuning parameters in (4) and (6). This resulted in values for
$(\mu_{i,0},\kappa_{i,0},\alpha_{i,0},\beta_{i,0})^{\top}$ that are listed in
Table 1. These tuning parameter values were used for both the CCP-PPM and that
of Loschi method.

Series | $\mu_{i,0}$ | $\kappa_{i,0}$ | $\alpha_{i,0}$ | $\beta_{i,0}$
---|---|---|---|---
USA | 0.009 | 188.924 | 2.212 | 1.233
Mexico | -0.014 | 9.075 | 1.596 | 0.574
Argentina | 0.009 | 10.186 | 1.328 | 0.402
Chile | -0.022 | 2.533 | 1.864 | 0.656
Brazil | 0.022 | 11.010 | 1.349 | 0.423
Table 1: Values of $(\mu_{i,0},\kappa_{i,0},\alpha_{i,0},\beta_{i,0})^{\top}$
for each country’s series. The values are the result of applying the empirical
procedure described in 3.2 to each country’s returns.

To specify values for $\bm{\mu}_{0}$ and $\bm{\Sigma_{0}}$, we first set
$m^{0}_{1}=\cdots=m^{0}_{5}=n^{-1}$ (in the application $n=1309$) and used a
compound symmetry matrix for $\mathbf{S}_{0}$ with variance
$\sigma^{2}_{0}=m^{0}_{1}(1-m^{0}_{1})/n=n^{-3}(n-1)$ and correlation $0.5$.
This resulted in $\bm{\mu}_{0}=(-7.1762)\bm{1}_{5}$ and $\bm{\Sigma}_{0}$
being a compound symmetric matrix with variance $0.334$ and correlation $0.5$.
We set $\nu_{0}=3$. For the Loschi method $p_{i}\sim\mbox{Beta}(a,b)$ and we
set $a=1304.5$ and $b=681209.9$. These values were selected based on setting
the mean number of clusters a priori to 3.5 with a variance of 2.5. Both
methods were fit by collecting 1000 MCMC samples after discarding the first
10,000 as burn-in and thinning by 5 (i.e., 15,000 total samples were
collected). Both the CCP-PPM and Loschi’s method were fit using the ccp_ppm
and icp_ppm functions that are available in the ppmSuite-package that can be
found on CRAN.

There are two approaches that could used to estimate change points. The first
classifies points as change points if their posterior probability of being a
change point is greater than some pre-specified value. The second classifies
points as change points based on a partition estimate. We report both as both
require input from the user (pre-specified probability cut-off for the first
approach and a loss function in the second approach).

We first explore the a posteriori dependence between partitions from the five
markets. To do this, at each MCMC iteration we computed the ARI for all
possible pairs of partitions (which is 10 in this application). The CCP-PPM
produced slightly more similar partitions across countries than the Loschi
method. The overall average pairwise ARI for CCP-PPM turned out to be 0.51
compared to 0.48 from the Loschi method.

Next we explore the posterior change point probabilities which are displayed
in the first column of Figure 5. The black points correspond to the CCP-PPM
and the red to Loschi. For both, change point probabilities were estimated
using the posterior means of $c_{i,t}$. It seems that there is a general
agreement between the two methods regarding the location of potential change
points. However, the CCP-PPM seems to produce probabilities that are closer to
one for these points compared to the Loschi method. In fact, the Loschi method
never records a change point probability greater than 0.75. Similarly, both
methods agree on the general location of points that have a small chance of
being a change point, although the CCP-PPM seems to push these probabilities
closer to zero compared to the Loschi method.

Figure 5: The left plot display posterior probability of each point being a
change point with black points corresponding to the CCP-PPM method and red
points to the Loschi method. The middle plot displays the partition estimate
of each series for the CCP-PPM and the right plot that for Loschi. Both sets
of partitions were estimated using the salso package in R.

Figure 5, second and third columns, shows the partition estimates under the
CCP-PPM (second column) and Loschi method (third column). Partition estimates
were obtained using the salso (Dahl et al., 2021b) R package and the
generalization of the Variation of Information loss function (Meilǎ, 2007).
Since in our case it seems natural to penalize change point false positives
more than false negatives, we set the false positive penalty parameter of the
salso function to $a=25$ (see Dahl et al., 2021a, for more details). We note
briefly that setting $a=25$ was driven primarily by Loschi’s method. If
$a>25$, then Loschi’s method tended to smooth over some change points and for
$a<25$ it tended to produce more change points than what would be desired.
However, the CCP-PPM was reasonably robust to $a$’s value. This is a
consequence of the change point probabilities from the the Loschi method being
more central (i.e., closer to 0.5) than those from the CCP-PPM.

Apparent differences between the estimated partitions exist, and they
illustrate how the CCP-PPM takes into account the dependency between the index
series or the contagion phenomenon. For example, the CCP-PPM method identifies
a shared partition for Brazil, Chile, Mexico, and the USA at the end of 1997.
It is important to stress that in July 1997, the Thai government ran out of
foreign currency, forcing it to float the Thai baht which is a factor in
starting the 1997 Asian financial crisis, or Asian Flu. This crisis spread
internationally, affecting some Asian stock markets such as Indonesia, South
Korea, Hong Kong, Laos, Malaysia, Philippines, Brunei, mainland China,
Singapore, Taiwan, and Vietnam. According to Harrigan (2000), the Asian
crises’ overall effect on the United States were small. However, as mentioned
by Stallings (1998), the Asian Flu hit the Latin American markets in October
1997, when bone spreads widened abruptly implying more risk. In the
Argentinean market, our approach identified a first partition at the end of
1996. Note that, at the end of 1995, the (real) gross domestic product (GDP)
in this country fell by 2.5 percent. However, by the end of 1996, it rebounded
by 5.5 percent (International Monetary Found, 2003), possibly affecting the
performance of the MERVAL index.

Another cluster or change point our method identifies is related to the
Russian crisis or Russian Cold in August 1998. This crisis started when the
Russian government and the Russian Central Bank devalued the ruble and
defaulted on its debt. It is important to stress that although most countries
experienced changes in their stock market returns series at the end of 1998,
the Argentinean market experienced a change in mid-1998. Moreover, the
IBOVESPA index (Brazil) experienced a change at the beginning of 1999, just
after the Russian crisis. This crisis was known as the Samba effect and was
produced when the Minas Gerais State Governor, Itamar Franco, stopped paying
Minas Gerais debt to other states, generating unleashing capital flights. Note
that a small cluster is detected by the CCP-PPM in the indexes of Argentina,
Mexico, and USA towards the end of 1998, but something the Loschi method
misses. These clusters provide some evidence that the contagion phenomena from
a mature market to emerging ones is being captured by the CCP-PPM, which
includes dependency between partitions.

Finally, the CCP-PPM also detected a cluster after the year 2000. The
corresponding change point can be explained by the dot-com bubble, caused by
excessive speculation of some internet companies in the late 1990s. On January
14, 2000, the Dow Jones Industrial Average reached its dot-com bubble peak.
This partition is observed in the Argentinean, Chilean, Mexican, and USA
markets. Possibly, the dot-com bubble may have affected other economies over
the Latin American region, evidencing some contagion effect. In this case,
Loschi method did not detect a change point at or near the above-mentioned
date.

Figure 6: Change point posterior probabilities. Gray points correspond to
locations whose posterior probability of being a change point was less than
0.4. The left column displays results under the CCP-PPM and the right column
under the Loschi method.

As mentioned, the estimated partitions in Figure 5 depend on the value $a=25$.
To show a more complete picture of both methods performance, we provide Figure
6. In this figure, all points with a posterior probability of being a change
point less than 0.4 are colored gray. The left column corresponds to the CCP-
PPM fit while the right the Loschi method. The CCP-PPM fit has more power in
detecting change points compared to Loschi method, without inflating the
false-positive rate. The points highlighted by the CCP-PPM fit are at least
plausibly change points, and those associated with more pronounced volatility
generally have a larger probability of being a change point, which is a
desirable characteristic.

##  6 Discussion

In this paper we developed a new change point detection model for $L$ time
series with $n$ observations in an arbitrary space, which undergo sudden
changes in their distributional parameters. By making dependent the vector of
latent change point probabilities $\bm{p}_{t}=(p_{1,t},\ldots,p_{L,t})^{\top}$
at a specific time $t\in\\{1,\ldots,n-1\\}$, the corresponding $L$ partitions
with contiguous clusters are encouraged to be correlated. We provide some
theoretical results that help to better understand the main features of our
model, a useful procedure to guide the specification of all parameters that
are involved in, and simple pseudo-code to perform posterior inference via
MCMC methods. Through a small simulation study, we compared the ability of our
model with other compelling approaches to detect highly dependent changes
under different scenarios, showing an improvement in detecting change points.
Additionally, we applied our method to the returns of emerging Latin American
and US markets, obtaining exciting results about possible contagion effect
between the economies of these countries based on the dependence between
change point locations.

In terms of extending the proposed model with the aim of making it more
flexible, several directions can be pursued. For instance, the assumption that
the vectors of change point probabilities
$\bm{p}_{t}=(p_{1,t},\ldots,p_{L,t})^{\top}$ are independent and identically
distributed through time $t$ can be relaxed. One possible approach would be to
model $\bm{p}_{t}$ with a stationary process in $(0,1)^{L}$. Another
interesting direction would be to incorporate time-dependent covariates in the
marginal likelihood function to describe abrupt changes in a regression curve.
A similar situation, but more complex, is to incorporate covariates in the
distribution for contiguous partitions. Finally, the computational cost
involved in the MCMC algorithm for posterior inference increases rapidly as
the length and number of time series grow. It would be appealing to develop
strategies that mitigate the so-called “curse of dimensionality”. These are
all topics of future research.

Acknowledgements

José J. Quinlan gratefully recognizes the financial support provided by the
Agencia Nacional de Investigación y Desarrollo de Chile (ANID) through
Fondecyt Grant 3190324.

##  7 Appendix

##  Appendix A Proofs of Propositions

Proofs of the propositions depend on a well known result of the multivariate
$t$. For the sake of completeness, we include it here.

######  Result 1 (Ding (2016)).

Let $A$ be a non-empty subset of $\\{1,\ldots,L\\}$. If
$\bm{z}=(z_{1},\ldots,z_{L})^{\top}\sim
t_{L}(\nu_{0},\bm{\mu}_{0},\bm{\Sigma}_{0})$, then $\bm{z}_{A}=(z_{l}:l\in
A)^{\top}\sim t_{|A|}(\nu_{0},\bm{\mu}_{A,0},\bm{\Sigma}_{A,0})$, where $|A|$
is the cardinality of $A$, $\bm{\mu}_{A,0}=(\mu^{0}_{l}:l\in A)^{\top}$ and
$\bm{\Sigma}_{A,0}=(\sigma^{0}_{l_{1},l_{2}}:l_{1},l_{2}\in A)$. When
$A=\\{a\\}$, set $z_{A}=z_{a}$, $\bm{\mu}_{A,0}=\mu^{0}_{a}$ and
$\bm{\Sigma}_{A,0}=\sigma^{0}_{a,a}$.

### Proof of Proposition 1.

A direct consequence of (6) and Result 1 is that

$\mathrm{logit}(p_{i,1}),\ldots,\mathrm{logit}(p_{i,n-1})\stackrel{{\scriptstyle\mathrm{iid.}}}{{\sim}}t_{1}(\nu_{0},\mu^{0}_{i},\sigma^{0}_{i,i}).$

Invoking the Change of Variables theorem, the marginal distribution for
$\bm{c}_{i}$ is given by

$\displaystyle\pi(\bm{c}_{i})$
$\displaystyle=\mathrm{E}\left[\prod_{t=1}^{n-1}p_{i,t}^{c_{i,t}}(1-p_{i,t})^{1-c_{i,t}}\right]$
$\displaystyle=\prod_{t=1}^{n-1}\mathrm{E}(p_{i,t})^{c_{i,t}}\\{1-\mathrm{E}(p_{i,t})\\}^{1-c_{i,t}}$
$\displaystyle=\prod_{t=1}^{n-1}\phi_{i}^{c_{i,t}}(1-\phi_{i})^{1-c_{i,t}},$

where

$\phi_{i}=\int_{\mathbb{R}}\left\\{\frac{\exp(z)}{1+\exp(z)}\right\\}t_{1}(z\mid\nu_{0},\mu^{0}_{i},\sigma^{0}_{i,i}).$

Keeping in mind that $k_{i}=1+\sum_{t=1}^{n-1}c_{i,t}$,

$\displaystyle\pi(\bm{c}_{i})$
$\displaystyle=\prod_{t=1}^{n-1}\phi_{i}^{c_{i,t}}(1-\phi_{i})^{1-c_{i,t}}$
$\displaystyle=\phi_{i}^{k_{i}-1}(1-\phi_{i})^{n-k_{i}}.$

From the last expression, the result follows.

### Proof of Proposition 2.

From the proof of Proposition 1, we saw that $c_{i,1},\ldots,c_{i,n-1}$ are
independent $\mathrm{Bernoulli}(\phi_{i})$ random variables. From this
observation, it follows that

$k_{i}-1=\sum_{t=1}^{n-1}c_{i,t}\sim\mathrm{Binomial}(n-1,\phi_{i}).$

As for the distribution of $(k_{i}-1,k_{s}-1)^{\top}$, given
$\bm{p}_{i,s}=(p_{i,1},p_{s,1},\ldots,p_{i,n-1},p_{s,n-1})^{\top}$ the random
variables $(k_{i}-1)$ and $(k_{s}-1)$ are conditionally independent Poisson’s
Binomial with corresponding parameters
$\bm{\tilde{p}}_{i}=(p_{i,1},\ldots,p_{i,n-1})^{\top}$ and
$\bm{\tilde{p}}_{s}=(p_{s,1},\ldots,p_{s,n-1})^{\top}$. The mixture
distribution follows since

$\displaystyle\mathbb{P}(k_{i}-1=a,k_{s}-1=b)$
$\displaystyle=\mathrm{E}\\{\mathbb{P}(k_{i}-1=a,k_{s}-1=b\mid\bm{p}_{i,s})\\}$
$\displaystyle=\mathrm{E}\\{\mathbb{P}(k_{i}-1=a\mid\bm{\tilde{p}}_{i})\mathbb{P}(k_{s}-1=b\mid\bm{\tilde{p}}_{s})\\},$

for all $a,b\in\\{0,\ldots,n-1\\}$. Now, by the Law of Total Covariance

$\displaystyle\mathrm{Cov}(k_{i}-1,k_{s}-1)$
$\displaystyle=\mathrm{Cov}\left(\sum_{t_{1}=1}^{n-1}c_{i,t_{1}},\sum_{t_{2}=1}^{n-1}c_{s,t_{2}}\right)$
$\displaystyle=\sum_{t_{1}=1}^{n-1}\sum_{t_{2}=1}^{n-1}\mathrm{Cov}(c_{i,t_{1}},c_{s,t_{2}})$
$\displaystyle=\sum_{t_{1}=1}^{n-1}\sum_{t_{2}=1}^{n-1}\mathrm{Cov}(p_{i,t_{1}},p_{s,t_{2}}).$

Due to (6) and Result 1 with its respective notation,

$(\mathrm{logit}(p_{i,1}),\mathrm{logit}(p_{s,1}))^{\top},\ldots,(\mathrm{logit}(p_{i,n-1}),\mathrm{logit}(p_{s,n-1}))^{\top}\stackrel{{\scriptstyle\mathrm{iid.}}}{{\sim}}t_{2}(\nu_{0},\bm{\mu}_{A,0},\bm{\Sigma}_{A,0}),$

where $A=\\{i,s\\}$. By the Change of Variables theorem and the previous
observation, we have that

$\mathrm{Cov}(p_{i,t_{1}},p_{s,t_{2}})=\mathrm{E}(p_{i,t_{1}}p_{s,t_{2}})-\mathrm{E}(p_{i,t_{1}})\mathrm{E}(p_{s,t_{2}})=\left\\{\begin{array}[]{lc}0,&\mathrm{if}\quad
t_{1}\neq t_{2}\\\ \varphi_{i,s}-\phi_{i}\phi_{s},&\mathrm{if}\quad
t_{1}=t_{2},\end{array}\right.$

where

$\varphi_{i,s}=\int_{\mathbb{R}^{2}}\left[\frac{\exp(z_{i})\exp(z_{s})}{\\{1+\exp(z_{i})\\}\\{1+\exp(z_{s})\\}}\right]t_{2}(\bm{z}_{A}\mid\nu_{0},\bm{\mu}_{A,0},\bm{\Sigma}_{A,0})\mathrm{d}\bm{z}_{A}.$

Finally,

$\mathrm{Cov}(k_{i}-1,k_{s}-1)=\sum_{t=1}^{n-1}\mathrm{Cov}(p_{i,t},p_{s,t})=(n-1)(\varphi_{i,s}-\phi_{i}\phi_{s}).$

### Proof of Proposition 3.

Consider two disjoint non-empty subsets $A,B$ of $\\{1,\ldots,L\\}$ and set
$\bm{c}_{D,t}=(c_{i,t}:i\in D)^{\top}$ for any non-empty subset $D$ of
$\\{1,\ldots,L\\}$. Then,

$\pi(\bm{c}_{A,t}\mid\bm{c}_{B,t})=\frac{\pi(\bm{c}_{A\cup
B,t})}{\pi(\bm{c}_{B,t})}.$

Using the Change of Variables theorem,

$\pi(\bm{c}_{B,t})=\int_{\mathbb{R}^{L}}\left\\{\prod_{i\in
B}\frac{\exp(c_{i,t}z_{i})}{1+\exp(z_{i})}\right\\}t_{L}(\bm{z}\mid\nu_{0},\bm{\mu}_{0},\bm{\Sigma}_{0})\mathrm{d}\bm{z},$

where $\bm{z}=(z_{1},\ldots,z_{L})^{\top}$. Using Result 1 and the notation
therein, it follows that

$\pi(\bm{c}_{B,t})=\int_{\mathbb{R}^{|B|}}\left\\{\prod_{i\in
B}\frac{\exp(c_{i,t}z_{i})}{1+\exp(z_{i})}\right\\}t_{|B|}(\bm{z}_{B}\mid\nu_{0},\bm{\mu}_{B,0},\bm{\Sigma}_{B,0})\mathrm{d}\bm{z}_{B}.$

Similar expression applies for $\pi(\bm{c}_{A\cup B,t})$. Therefore,

$\pi(\bm{c}_{A,t}\mid\bm{c}_{B,t})=\frac{\int_{\mathbb{R}^{|A|+|B|}}\left\\{\prod_{i\in
A\cup
B}\frac{\exp(c_{i,t}z_{i})}{1+\exp(z_{i})}\right\\}t_{|A|+|B|}(\bm{z}_{A\cup
B}\mid\nu_{0},\bm{\mu}_{A\cup B,0},\bm{\Sigma}_{A\cup
B,0})\mathrm{d}\bm{z}_{A\cup B}}{\int_{\mathbb{R}^{|B|}}\left\\{\prod_{i\in
B}\frac{\exp(c_{i,t}z_{i})}{1+\exp(z_{i})}\right\\}t_{|B|}(\bm{z}_{B}\mid\nu_{0},\bm{\mu}_{B,0},\bm{\Sigma}_{B,0})\mathrm{d}\bm{z}_{B}}.$

As for the conditional probability $\mathbb{P}(c_{i,t}=1\mid c_{s,t}=1)$, take
$A=\\{i\\}$ and $B=\\{s\\}$.

## References

  * Arellano-Valle et al. (2013) Arellano-Valle, R., Castro, L. M., and Loschi, R. (2013), “Change point detection in the skew-normal model parameters,” Communications in Statistics. Theory and Methods, 42, 603–618.
  * Barry and Hartigan (1992) Barry, D. and Hartigan, J. A. (1992), “Product partition models for change point problems,” Ann. Statist., 20, 260–279.
  * Cheon and Kim (2010) Cheon, S. and Kim, J. (2010), “Multiple Change-Point Detection of Multivariate Mean Vectors with the Bayesian Approach,” Comput. Stat. Data Anal., 54, 406––415.
  * Dahl et al. (2021a) Dahl, D. B., Johnson, D. J., and Mueller, P. (2021a), “Search Algorithms and Loss Functions for Bayesian Clustering,” arXiv:2105.04451v1.
  * Dahl et al. (2021b) Dahl, D. B., Johnson, D. J., and Müller, P. (2021b), salso: Search Algorithms and Loss Functions for Bayesian Clustering, r package version 0.2.23.
  * de P. Filleti et al. (2008) de P. Filleti, J., Hotta, L. K., and Zevallos, M. (2008), “Analysis of contagion in emerging markets,” Journal of Data Science, 6, 601–626.
  * Ding (2016) Ding, P. (2016), “On the Conditional Distribution of the Multivariate t Distribution,” The American Statistician, 70, 293–295.
  * Erdman and Emerson (2007) Erdman, C. and Emerson, J. W. (2007), “bcp: An R Package for Performing a Bayesian Analysis of Change Point Problems,” Journal of Statistical Software, 23, 1–13.
  * Fan and Mackey (2017) Fan, Z. and Mackey, L. (2017), “Empirical Bayesian Analysis of Simultaneous Changepoints in Multiple Data Sequences,” Annals of Applied Statistics, 11, 2200–2221.
  * García and Gutiérrez-Peña (2019) García, E. C. and Gutiérrez-Peña, E. (2019), “Nonparametric product partition models for multiple change-points analysis,” Communications in Statistics - Simulation and Computation, 48, 1922–1947.
  * Gupta et al. (2021) Gupta, S. K., Gupta, N., and Singh, V. P. (2021), “Variable-Sized Cluster Analysis for 3D Pattern Characterization of Trends in Precipitation and Change-Point Detection,” Journal of Hydrologic Engineering, 26, 04020056\.
  * Harlé et al. (2016) Harlé, F., Chatelain, F., Gouy-Pailler, C., and Achard, S. (2016), “Bayesian Model for Multiple Change-Points Detection in Multivariate Time Series,” IEEE Transactions on Signal Processing, 64, 4351–4362.
  * Harrigan (2000) Harrigan, J. (2000), “The impact of the Asia crisis on U.S industry: An almost-free lunch?” Federal Reserve Bank of New York, Economic Policy Review.
  * International Monetary Found (2003) International Monetary Found (2003), “El papel del FMI en la Argentina, 1991-2002,” .
  * James et al. (2019) James, N. A., Zhang, W., and Matteson, D. S. (2019), “ecp: An R Package for Nonparametric Multiple Change Point Analysis of Multivariate Data. R package version 3.1.2,” .
  * Jin et al. (2021) Jin, H., Yin, G., Yuan, B., and Jiang, F. (2021), “Bayesian Hierarchical Model for Change Point Detection in Multivariate Sequences,” Technometrics, 0, 1–30.
  * Jones et al. (2021) Jones, C., Clayton, S., Ribalet, F., Armbrust, E. V., and Harchaoui, Z. (2021), “A kernel-based change detection method to map shifts in phytoplankton communities measured by flow cytometry,” Methods in Ecology and Evolution, 12, 1687–1698.
  * Killick et al. (2012) Killick, R., Fearnhead, P., and Eckley, I. A. (2012), “Optimal Detection of Changepoints With a Linear Computational Cost,” Journal of the American Statistical Association, 107, 1590–1598.
  * Loschi and Cruz (2002) Loschi, R. and Cruz, F. (2002), “Analysis of the influence of some prior specifications in the identification of change points via product partition model,” Computational Statistics & Data Analysis, 39, 477–501.
  * Loschi and Cruz (2005) — (2005), “Extension to the product partition model: computing the probability of a change,” Computational Statistics & Data Analysis, 48, 255–268.
  * Loschi et al. (2005) Loschi, R., Cruz, F., and Arellano-Valle, R. (2005), “Multiple change point analysis for the regular exponential family using the product partition model,” Journal of Data Science, 3, 305–330.
  * Loschi et al. (2003) Loschi, R., Cruz, F., Iglesias, P., and Arellano-Valle, R. (2003), “A Gibbs sampling scheme to product partition model: An application to change-point problems,” Computers & Operations Research, 30, 463–482.
  * Loschi et al. (2010) Loschi, R., Pontel, J., and Cruz, F. (2010), “Multiple change-point analysis for linear regression models,” Chilean Journal of Statistics, 1, 93–112.
  * Lowell et al. (1998) Lowell, J. F., Neu, C. R., and Tong, D. (1998), Financial Crises and Contagion in Emerging Market Countries, Santa Monica, CA: RAND Corporation.
  * Martínez and Mena (2014) Martínez, A. F. and Mena, R. H. (2014), “On a Nonparametric Change Point Detection Model in Markovian Regimes,” Bayesian Analysis, 9, 823–858.
  * Matteson and James (2014) Matteson, D. and James, N. (2014), “A Nonparametric Approach for Multiple Change Point Analysis of Multivariate Data,” Journal of the American Statistical Association, 109, 334–345.
  * Meilǎ (2007) Meilǎ, M. (2007), “Comparing clusterings-an information based distance,” Journal of Multivariate Analysis, 98, 873–895.
  * Nyamundanda et al. (2015) Nyamundanda, G., Hegarty, A., and Hayes, K. (2015), “Product partition latent variable model for multiple change-point detection in multivariate data,” Journal of Applied Statistics, 42, 2321–2334.
  * Padilla et al. (2021) Padilla, O. H. M., Yu, Y., Wang, D., and Rinaldo, A. (2021), “Optimal nonparametric change point analysis,” Electronic Journal of Statistics, 15, 1154–1201.
  * Page and Quinlan (2022) Page, G. L. and Quinlan, J. J. (2022), ppmSuite: A Collection of Models that Employ a Product Parition Prior Distribution on Partitions, r package version 0.2.1.
  * Page et al. (2021) Page, G. L., Quintana, F. A., and Dahl, D. B. (2021), “Dependent Modeling of Temporal Sequences of Random Partitions,” Journal of Computational and Graphical Statistics, 0, 1–29.
  * Pedroso et al. (2021) Pedroso, R. C., Loschi, R. H., and Quintana, F. A. (2021), “Multipartition model for multiple change point identification,” arXiv:2107.11456v2.
  * R Core Team (2021) R Core Team (2021), R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria.
  * Rand (1971) Rand, W. M. (1971), “Objective Criteria for the Evaluation of Clustering Methods,” Journal of the American Statistical Association, 66, 846–850.
  * Stallings (1998) Stallings, B. (1998), “Impact of th Asian crisis on Latin America,” United Nations Economic Commission for Latin America and the Caribbean.
  * Truong et al. (2020) Truong, C., Oudre, L., and Vayatis, N. (2020), “Selective review of offline change point detection methods,” Signal Processing, 167, 107299\.
  * Tveten et al. (2021) Tveten, M., Eckley, I. A., and Fearnhead, P. (2021), “Scalable changepoint and anomaly detection in cross-correlated data with an application to condition monitoring,” arXiv:2010.06937v2.
  * Valdés (2000) Valdés, R. (2000), “Emerging markets contagion: evidence and theory,” Available at SSRN: https://ssrn.com/abstract=69093 or http://dx.doi.org/10.2139/ssrn.69093.
  * Wang and Emerson (2015) Wang, X. and Emerson, J. W. (2015), “Bayesian Change Point Analysis of Linear Models on Graphs,” arXiv:1509.00817.
  * Wang (1993) Wang, Y. H. (1993), “On the Number of Successes in Independent Trials,” Statistica Sinica, 3, 295–312.
  * Wood et al. (2021) Wood, K., Roberts, S., and Zohren, S. (2021), “Slow Momentum with Fast Reversion: A Trading Strategy Using Deep Learning and Changepoint Detection,” The Journal of Financial Data Science, jfds.2021.1.081.
  * Yao (1984) Yao, Y.-C. (1984), “Estimation of a noisy discrete-time step function: Bayes and empirical Bayes approaches,” Ann. Statist., 12, 1434–1447.
  * Zanini et al. (2019) Zanini, C. T. P., Müller, P., Ji, Y., and Quintana, F. A. (2019), “A Bayesian Random Partition Model for Sequential Refinement and Coagulation,” Biometrics, 75, 988–999.