# Design a Delicious Lunchbox in Style

Yutong Zhou
Ritsumeikan University, Shiga, Japan
<EMAIL_ADDRESS>
###### Abstract

We propose a cyclic generative adversarial network with spatial-wise and
channel-wise attention modules for text-to-image synthesis [8]. To accurately
depict and design scenes with multiple occluded objects, we design a pre-
trained ordering recovery model and a generative adversarial network to
predict layout and composite novel box lunch presentations [9]. In the
experiments, we devise the Bento800 dataset to evaluate the performance of the
text-to-image synthesis model and the layout generation & image composition
model. This paper is the continuation of our previous paper works [8] and [9].
We also present additional experiments and qualitative performance comparisons
to verify the effectiveness of our proposed method. Bento800 dataset is
available at https://github.com/Yutong-Zhou-cv/Bento800_Dataset.

##  1 Introduction

With the increasing need for portable meals as people had to leave their homes
for hours, the bento box has become a popular and celebrated food culture that
emphasizes flavor, visual aesthetics, nutritional balance, and convenience.
Despite advances in robotics applications in food production, creating
visually appealing box lunch presentations and organizing multiple food items
in an orderly manner still remain a challenge for robotics. This raises a
question: What kind of composite image can provide comprehensive robotics
guidance, standard placement compliance, and visually appealing presentation?

To resolve the above issue, (1) We propose a cyclic generative adversarial
network [8] for text-to-image generation and image captioning. (2) We
introduce Bento800 [9], the first manually annotated synthetic box lunch
dataset. (3) We propose an aesthetic box lunch design model [9] with pre-
trained placement ordering recovery and a generative adversarial network (GAN)
[2] for layout generation & ingredients composition. This paper effectively
combines our prior works and additionally expands upon the following aspects
(Details are shown in Sec. 2.3 and Sec. 3): (4) We rephrase text descriptions
in Bento800 dataset to increase text diversity for text-to-image generation.
(5) We conduct experiments to demonstrate the effectiveness of our method in
synthesizing and designing novel lunchbox images.

##  2 Method

###  2.1 Text-to-Image Synthesis

Figure 1: An overview of our proposed approach in [8].

Based on our previous dual attention architecture [7], we propose a cyclic GAN
(as in Fig. 1) and adopt a pre-trained BERT base model [1] for text embedding.
For image captioning, we employ the Inception-v3 model [6] as a CNN image
encoder for extracting image features and the BERT model to generate
contextualized word vectors for each caption. The full objective function as
below:

${L}_{G}=\sum_{m=1}^{M}({L}_{G_{m}}+{L}_{id})+{L}_{cycle}$ (1)
${L}_{D}=\sum_{m=1}^{M}{L}_{D_{m}}$ (2)

where ${L}_{G_{m}}$ and ${L}_{D_{m}}$ denote the adversarial loss at mth stage
(m=0, 1, 2, 3). The adversarial loss consists of two components: unconditional
loss, which evaluates the realism of the generated image, and conditional
loss, which measures the compatibility between the generated image and the
input text description. ${L}_{id}$ denotes the proposed identity-preserving
loss which quantifies the dissimilarity between the generated image and the
ground-truth image. ${L}_{cycle}$ denotes the cycle consistency loss between
predicted captions and target captions.

###  2.2 Layout Generation and Image Composition

Text-guided generation models often have imprecise specifications for
producing the intended target image, necessitating the use of additional
controls for generating images. To achieve a precise depiction of scenes with
multiple objects, we propose the following method, as shown in Fig. 2. The
pre-processing model in Fig. 2(a) takes synthetic or real images as input and
builds a placement ordering list. During training in Fig. 2(a), the layout
generation model predicts layout bounding boxes while the ingredients
composition model transforms food items to their corresponding transposed
targets. During the testing in Fig. 2(b), the weights of the layout generation
model, Spatial Transformer Network (STN) [3] model and ingredients composition
model are frozen. The training objective function is computed as:

$\centering{L}={L}_{layout}(G,D)+{L}_{image}(G,D)+{L}_{STN}\vspace{-1mm}\@add@centering$
(3)

Conditional GAN loss functions (${L}_{layout}$ and ${L}_{image}$) are applied
for Layout Generation and Ingredient composition network. We also adopt a
differentiable grid generator and bilinear sampling to transform food items to
proper positions and compute the ${L}_{STN}$ (L1 loss) between input single-
item food and transformed single-item food.

(a) Pre-processing & Training Phase

(b) Testing Phase

Figure 2: An overview of our proposed approach in [9].

###  2.3 Extend Information: Paraphrasing

We propose a novel dataset Bento800 [9], which comprises 800 lunchbox images
with three different types of food presentations: (1) Place fried chicken on
rice; (2) Place salt-grilled salmon and tamagoyaki on rice; (3) Place
croquette and fried shrimp on rice, fried shrimp is on croquette. To produce
parallel and diverse data, we use OpenAI chatGPT111https://chat.openai.com/ to
rephrase the above sentences with the template “`Rephrase: [Input] with 25
examples`” and “`Rephrase: [Input]`”+“`Please try again`”$\times 24$ and also
perform data cleaning. Then, we randomly select 8 sentences from the rephrased
results for each lunchbox image and merge them with the original sentences to
expand our dataset for text-to-image generation.

##  3 Experiments

Figure 3: Qualitative comparison: Ours [8] (left), Composable-Diffusion [4]
(middle) and Stable Diffusion [5] (right). Figure 4: Our synthesized box lunch
presentations [9] (“Ours”) and human-designed ground-truth (“GT”). Figure 5:
Diverse box lunch presentations [9].

Figure 3 illustrates a qualitative comparison between box lunch presentations
with state-of-the-art methods for text-to-image. However, they can only
synthesize poor image quality results [8] or images corresponding to part of
the text prompt [4, 5]. Figure 4 demonstrates more accurate and high-quality
box lunch presentations organized by our proposed method [9] and human-
designed. Three groups of results show different types of food presentations,
which are also present in Bento800. Figure 5 shows various lunchboxes from
single food items by randomly sampling noise. Experimental results show that
our model produces diverse lunchboxes that conform to popular aesthetic
preferences.

##  4 Conclusion

This paper presents a cyclic GAN for text-to-image synthesis and a layout
generation & food composition network for innovative lunchbox presentation
design. We also introduce Bento800, the first manually annotated dataset with
extensive annotations. For future work, we plan to create a metric that aligns
with human aesthetic judgment. This exploration may encourage the potential of
more varied applications for bento box presentation design in robotics.

## References

  * [1] Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  Bert: Pre-training of deep bidirectional transformers for language understanding.  arXiv preprint arXiv:1810.04805, 2018.
  * [2] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial nets.  Advances in neural information processing systems, 27, 2014.
  * [3] Max Jaderberg, Karen Simonyan, Andrew Zisserman, et al.  Spatial transformer networks.  Advances in neural information processing systems, 28, 2015.
  * [4] Nan Liu, Shuang Li, Yilun Du, Antonio Torralba, and Joshua B Tenenbaum.  Compositional visual generation with composable diffusion models.  In Computer Vision–ECCV 2022: 17th European Conference, Tel Aviv, Israel, October 23–27, 2022, Proceedings, Part XVII, pages 423–439. Springer, 2022.
  * [5] Robin Rombach, Andreas Blattmann, Dominik Lorenz, Patrick Esser, and Björn Ommer.  High-resolution image synthesis with latent diffusion models.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 10684–10695, 2022.
  * [6] Christian Szegedy, Vincent Vanhoucke, Sergey Ioffe, Jon Shlens, and Zbigniew Wojna.  Rethinking the inception architecture for computer vision.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), June 2016.
  * [7] Yutong Zhou.  Generative adversarial network for text-to-face synthesis and manipulation.  In Proceedings of the 29th ACM International Conference on Multimedia, pages 2940–2944, 2021.
  * [8] Yutong Zhou and Nobutaka Shimada.  Generative adversarial network for text-to-face synthesis and manipulation with pretrained bert model.  In 2021 16th IEEE International Conference on Automatic Face and Gesture Recognition (FG 2021), pages 01–08. IEEE, 2021.
  * [9] Yutong Zhou and Nobutaka Shimada.  Able: Aesthetic box lunch editing.  In Proceedings of the 1st International Workshop on Multimedia for Cooking, Eating, and related APPlications, pages 53–56, 2022.


$\displaystyle=\widetilde{\textbf{G}}_{c}(\pi)-\widetilde{\textbf{G}}_{c}(\pi)\bigg{[}\widetilde{\textbf{G}}_{c}(\pi)+\widetilde{\textbf{G}}_{c}(1)-\widetilde{\textbf{G}}_{c}(\pi)\bigg{]}^{-1}\widetilde{\textbf{G}}_{c}(\pi)$
$\displaystyle=\widetilde{\textbf{G}}_{c}(\pi)-\widetilde{\textbf{G}}_{c}(\pi)\widetilde{\textbf{{G}}}_{c}(1)^{-1}\widetilde{\textbf{G}}_{c}(\pi)$

Thus, we show that the limiting distribution of the sup OLS-Wald statistic is
given by

$\displaystyle\widetilde{\mathcal{W}}^{OLS}(\pi)\equiv\underset{\pi\in[\pi_{1},\pi_{2}]}{\text{
sup }}\
\mathcal{W}_{T}^{OLS}(\pi)\Rightarrow\displaystyle\underset{\pi\in[\pi_{1},\pi_{2}]}{\text{
sup }}\ \
\bigg{\\{}\widetilde{\mathbf{N}}^{\prime}_{c}(\pi)\widetilde{\mathbf{M}}_{c}(\pi)^{-1}\widetilde{\mathbf{N}}_{c}(\pi)\bigg{\\}}$
(A.13)

with quantities $\widetilde{\mathbf{M}}_{c}(\pi)$,
$\widetilde{\mathbf{N}}_{c}(\pi)$,$\widetilde{\textbf{G}}_{c}(\pi)$ and
$\widetilde{H}_{c}(\pi)$ as defined by Proposition 1. ∎

Proof of Theorem 2.

###### Proof.

Let $Y$ denote the vector with all demeaned values of $y_{t}$ and $X$ be the
matrix collecting all demeaned values of $x_{t-1}$, that is,

$\displaystyle
Y=\left(y_{1}-\bar{y}_{T},y_{2}-\bar{y}_{T},...,y_{T}-\bar{y}_{T}\right)^{\prime}\text{and}\
X=\left(x_{0}-\bar{x}_{T-1},x_{1}-\bar{x}_{T-1},...,x_{T-1}-\bar{x}_{T-1}\right)^{\prime}.$

Similarly, we use $\mathcal{U}$ to denote the corresponding demeaned $u_{t}$
vector. Furthermore, for any $1\leq t\leq T$, we define $X_{t}$ to be a
$T\times p$ matrix, whose first $t$ rows are the same as $X$ while the rest
are all zeros. Moreover, let $Z=\left(z_{0},z_{1},...z_{T-1}\right)^{\prime}$
collect all the IVX instruments, and
$Z_{t}=\left(z_{0},...,z_{t},0,....,0\right)^{\prime}$ be the corresponding
time$-t$ truncated matrix. Given these notations, we express the original
predictive regression model as

$\displaystyle Y=X\beta_{2}+X_{t}\phi+\mathcal{U}$ (A.14)

where $\phi:=\beta_{2}-\beta_{1}$ measures the magnitude of structural break.
Moreover, we denote with $\phi_{t}$ to the corresponding estimator which
captures the break size associated with the sample partition at time $t$.
Therefore, given any particular $t$, testing for structural break in the
parameter vector $\beta$ is equivalent to testing the null hypothesis
$\phi_{t}=0$. Define with
$\mathbf{M}_{xz}=\mathbf{I}_{p}-X\left(Z^{\prime}X\right)^{-1}Z^{\prime}$,
which is idempotent and orthogonal to both $X$ and $Z$ and allows to rewrite
(A.14) in its canonical form303030Notice that the reparametrization of the
model to its canonical form allows to shift the coordinates which transforms
the model to its more general form within the exponential family.. Multiplying
$\mathbf{M}_{xz}$ on both sides of (A.14), we deduce that
$\mathbf{M}_{xz}Y=\mathbf{M}_{xz}X_{t}\phi_{t}+\mathbf{M}_{xz}\mathcal{U}$.
Now, using $\mathbf{M}_{xz}Z_{t}$ as the instrumental variables for
$\mathbf{M}_{xz}X_{t}$, we obtain an estimator for the parameter $\phi_{t}$
given by

$\displaystyle\widetilde{\phi}_{t}=\left(Z_{t}^{\prime}\mathbf{M}_{xz}X_{t}\right)^{-1}Z_{t}^{\prime}\mathbf{M}_{xz}Y$
(A.15)

Moreover, it holds that
$\widetilde{\phi}_{t}=\widetilde{\beta}_{2}^{IVX}-\widetilde{\beta}_{1}^{IVX}$.
Thus, the limiting distribution is given by

$\displaystyle\widetilde{\phi}_{t}-\phi_{t}$
$\displaystyle=\left(Z_{t}^{\prime}\mathbf{M}_{xz}X_{t}\right)^{-1}Z_{t}^{\prime}\mathbf{M}_{xz}\mathcal{U}$
$\displaystyle=\bigg{(}Z_{t}^{\prime}X_{t}-Z_{t}^{\prime}X\left(Z^{\prime}X\right)^{-1}Z^{\prime}X_{t}\bigg{)}^{-1}\bigg{(}Z_{t}^{\prime}\mathcal{U}_{y}-Z_{t}^{\prime}X_{t}\left(Z^{\prime}X\right)^{-1}Z^{\prime}\mathcal{U}\bigg{)}$
$\displaystyle=\left[\sum_{t=1}^{\lfloor
T\pi\rfloor}\widetilde{z}_{t}\left(x_{t}-\bar{x}_{T-1}\right)^{\prime}-\sum_{t=1}^{\lfloor
T\pi\rfloor}\widetilde{z}_{t}\left(x_{t}-\bar{x}_{T-1}\right)^{\prime}\left(\sum_{t=1}^{T}\widetilde{z}_{t}\left(x_{t}-\bar{x}_{T-1}\right)^{\prime}\right)^{-1}\sum_{t=1}^{\lfloor
T\pi\rfloor}\widetilde{z}_{t}\left(x_{t}-\bar{x}_{T-1}\right)^{\prime}\right]^{-1}$
$\displaystyle\times\left[\sum_{t=1}^{\lfloor
T\pi\rfloor}\widetilde{z}_{t}\left(u_{t}-\bar{u}_{T}\right)^{\prime}-\sum_{t=1}^{\lfloor
T\pi\rfloor}\widetilde{z}_{t}\left(x_{t}-\bar{x}_{T-1}\right)^{\prime}\left(\sum_{t=1}^{T}\widetilde{z}_{t}\left(x_{t}-\bar{x}_{T-1}\right)^{\prime}\right)^{-1}\sum_{t=1}^{T}\widetilde{z}_{t}\left(u_{t}-\bar{u}_{T}\right)\right]$

Applying the asymptotic results given by Lemma 1 we obtain

$\displaystyle
T^{\frac{1+\gamma_{x}\wedge\delta_{z}}{2}}\left(\widetilde{\phi}_{t}-\phi_{t}\right)\Rightarrow\bigg{[}\mathbf{\Psi}(\pi)-\mathbf{\Psi}(\pi)\mathbf{\Psi}(1)^{-1}\mathbf{\Psi}(\pi)^{\prime}\bigg{]}^{-1}\bigg{(}U(\pi)-\mathbf{\Psi}(\pi)\mathbf{\Psi}(1)^{-1}U(1)\bigg{)}$

Next, we focus on covariance estimators for
$\widetilde{\mathbf{Q}}_{1}(t)=\left(Z_{t}^{\prime}X_{t}\right)^{-1}\left(Z_{t}^{\prime}Z_{t}\right)\left(X_{t}^{\prime}Z_{t}\right)^{-1}$
and the corresponding one for $\widetilde{\mathbf{Q}}_{2}(t)$. By ignoring the
second-order degree bias correction we obtain:

$\displaystyle\widetilde{\mathbf{Q}}_{1}(t)=\left(\sum_{t=1}^{\lfloor
T\pi\rfloor}\widetilde{z}_{t}\left(x_{t}-\bar{x}_{T-1}\right)^{\prime}\right)^{-1}\left(\sum_{t=1}^{\lfloor
T\pi\rfloor}\widetilde{z}_{t}\widetilde{z}_{t}^{\prime}\right)\left(\sum_{t=1}^{\lfloor
T\pi\rfloor}\left(x_{t}-\bar{x}_{T-1}\right)\widetilde{z}_{t}^{\prime}\right)^{-1}$
(A.16)

which implies that

$\displaystyle
T^{1+\gamma_{x}\wedge\delta_{z}}\widetilde{\mathbf{Q}}_{1}(t)\Rightarrow\mathbf{\Psi}(\pi)^{-1}\left(\pi\sigma^{2}_{u}\widetilde{\mathbf{V}}\right)\mathbf{\Psi}(\pi)^{-1^{\prime}}$
(A.17)

Similarly,

$\displaystyle\widetilde{\mathbf{Q}}_{2}(t)=\left(\sum_{t=\lfloor
T\pi\rfloor+1}^{T}\widetilde{z}_{t}\left(x_{t}-\bar{x}_{T-1}\right)^{\prime}\right)^{-1}\left(\sum_{t=\lfloor
T\pi\rfloor+1}^{T}\widetilde{z}_{t}\widetilde{z}_{t}^{\prime}\right)\left(\sum_{t=\lfloor
T\pi\rfloor+1}^{T}\left(x_{t}-\bar{x}_{T-1}\right)\widetilde{z}_{t}^{\prime}\right)^{-1}$
(A.18)

which implies that

$\displaystyle
T^{1+\gamma_{x}\wedge\delta_{z}}\widetilde{\mathbf{Q}}_{2}(t)\Rightarrow\big{(}\mathbf{\Psi}(1)-\mathbf{\Psi}(\pi)\big{)}^{-1}\left((1-\pi)\sigma^{2}_{u}\widetilde{\mathbf{V}}\right)\big{(}\mathbf{\Psi}(1)-\mathbf{\Psi}(\pi)\big{)}^{-1^{\prime}}$
(A.19)

Combining all the above we obtain the following result for the
$\mathcal{W}_{b}(t)$ test statistic

$\displaystyle\mathcal{W}_{b}(t)$
$\displaystyle=\bigg{(}\widetilde{\beta}_{2}^{IVX}(t)-\widetilde{\beta}_{1}^{IVX}(t)\bigg{)}^{\prime}\left[\widetilde{\mathbf{Q}}_{1}(t)+\widetilde{\mathbf{Q}}_{2}(t)\right]^{-1}\bigg{(}\widetilde{\beta}_{2}^{IVX}(t)-\widetilde{\beta}_{1}^{IVX}(t)\bigg{)}$
$\displaystyle=\left\\{T^{\frac{1+\gamma_{x}\wedge\delta_{z}}{2}}\bigg{(}\widetilde{\beta}_{2}^{IVX}(t)-\widetilde{\beta}_{1}^{IVX}(t)\bigg{)}\right\\}^{\prime}\times\left[T^{1+\gamma_{x}\wedge\delta_{z}}\left(\widetilde{\mathbf{Q}}_{1}(t)+\widetilde{\mathbf{Q}}_{2}(t)\right)\right]^{-1}$
$\displaystyle\times\left\\{T^{\frac{1+\gamma_{x}\wedge\delta_{z}}{2}}\bigg{(}\widetilde{\beta}_{2}^{IVX}(t)-\widetilde{\beta}_{1}^{IVX}(t)\bigg{)}\right\\}$
(A.20)

which implies the following weakly convergence result

$\displaystyle\mathcal{W}^{IVX}_{\beta}(t)$
$\displaystyle\Rightarrow\bigg{(}U(\pi)-\mathbf{\Psi}(\pi)\mathbf{\Psi}(1)^{-1}U(1)\bigg{)}^{\prime}\bigg{[}\mathbf{\Psi}(\pi)-\mathbf{\Psi}(\pi)\mathbf{\Psi}(1)^{-1}\mathbf{\Psi}(\pi)^{\prime}\bigg{]}^{-1^{\prime}}$
$\displaystyle\times\bigg{[}\mathbf{\Psi}(\pi)^{-1}\left(\pi\sigma^{2}_{y}\widetilde{\mathbf{V}}\right)\mathbf{\Psi}(\pi)^{-1^{\prime}}+\big{(}\mathbf{\Psi}(1)-\mathbf{\Psi}(\pi)\big{)}^{-1}\left((1-\pi)\sigma^{2}_{y}\widetilde{\mathbf{V}}\right)\big{(}\mathbf{\Psi}(1)-\mathbf{\Psi}(\pi)\big{)}^{-1^{\prime}}\bigg{]}^{-1}$
$\displaystyle\times\bigg{[}\mathbf{\Psi}(\pi)-\mathbf{\Psi}(\pi)\mathbf{\Psi}(1)^{-1}\mathbf{\Psi}(\pi)^{\prime}\bigg{]}^{-1}\bigg{(}U(\pi)-\mathbf{\Psi}(\pi)\mathbf{\Psi}(1)^{-1}U(1)\bigg{)}$
(A.21)

To simplify the notation of expression (A), we denote with
$A=\mathbf{\Psi}(\pi)$, $C=\mathbf{\Psi}(1)$ and
$\Sigma=\sigma^{2}_{u}\widetilde{\mathbf{V}}$. Then, we have the following
equivalent form of the IVX-Wald statistic

$\displaystyle\mathcal{W}^{IVX}_{\beta}(t)\equiv{}$
$\displaystyle\bigg{(}U(\pi)-AC^{-1}U(1)\bigg{)}^{\prime}\bigg{(}A-AC^{-1}A^{\prime}\bigg{)}^{-1^{\prime}}$
(A.22) $\displaystyle\times\bigg{[}\pi A^{-1}\Sigma
A^{-1^{\prime}}+(1-\pi)\big{(}C-A\big{)}^{-1}\Sigma\big{(}C-A\big{)}^{-1^{\prime}}\bigg{]}^{-1}$
$\displaystyle\times\bigg{(}A-AC^{-1}A^{\prime}\bigg{)}^{-1}\bigg{(}U(\pi)-AC^{-1}U(1)\bigg{)}$

which can be written as below

$\displaystyle\mathcal{W}^{IVX}_{\beta}(t)\equiv{}$
$\displaystyle\bigg{(}U(\pi)-AC^{-1}U(1)\bigg{)}^{\prime}$ (A.23)
$\displaystyle\times\bigg{(}\pi\left[\bigg{(}A-AC^{-1}A^{\prime}\bigg{)}A^{-1}\Sigma
A^{-1^{\prime}}\bigg{(}A-AC^{-1}A^{\prime}\bigg{)}^{\prime}\right]$
$\displaystyle+(1-\pi)\left[\bigg{(}A-AC^{-1}A^{\prime}\bigg{)}\left(C-A\right)^{-1}\Sigma\left(C-A\right)^{-1^{\prime}}\bigg{(}A-AC^{-1}A^{\prime}\bigg{)}^{\prime}\right]\bigg{)}$
$\displaystyle\times\bigg{(}U(\pi)-AC^{-1}U(1)\bigg{)}$
$\displaystyle=\bigg{(}U(\pi)-AC^{-1}U(1)\bigg{)}^{\prime}$
$\displaystyle\times\bigg{(}\pi\big{(}I-AC^{-1}\big{)}\Sigma\big{(}I-AC^{-1}\big{)}^{\prime}+(1-\pi)\left(AC^{-1}\right)\Sigma\left(AC^{-1}\right)^{\prime}\bigg{)}^{-1}$
$\displaystyle\times\bigg{(}U(\pi)-AC^{-1}U(1)\bigg{)}.$

Notice that since $U(.)$ is known to be a Brownian motion with variance
$\Sigma$, the above expression can be further simplified as following

$\displaystyle\mathcal{W}^{IVX}_{\beta}(t)\equiv{}$
$\displaystyle\bigg{(}B(\pi)-AC^{-1}B(1)\bigg{)}^{\prime}$ (A.24)
$\displaystyle\times\bigg{(}\pi\big{(}I-AC^{-1}\big{)}\big{(}I-AC^{-1}\big{)}^{\prime}+(1-\pi)\left(AC^{-1}\right)\left(AC^{-1}\right)^{\prime}\bigg{)}^{-1}$
$\displaystyle\times\bigg{(}B(\pi)-AC^{-1}B(1)\bigg{)}.$

Substituting back to expression (A.24) the notation for
$A=\mathbf{\Psi}(\pi)$, $C=\mathbf{\Psi}(1)$ and
$\Sigma=\sigma^{2}_{u}\widetilde{\mathbf{V}}$, we obtain the expression below

$\displaystyle\mathcal{W}^{IVX}_{\beta}(t)\Rightarrow{}$
$\displaystyle\bigg{(}B(\pi)-\mathbf{\Psi}(\pi)\mathbf{\Psi}(1)^{-1}B(1)\bigg{)}^{\prime}$
(A.25)
$\displaystyle\times\bigg{(}\pi\big{(}\mathbf{I}_{p}-\mathbf{\Psi}(\pi)\mathbf{\Psi}(1)^{-1}\big{)}\big{(}\mathbf{I}_{p}-\mathbf{\Psi}(\pi)\mathbf{\Psi}(1)^{-1}\big{)}^{\prime}+(1-\pi)\left(\mathbf{\Psi}(\pi)\mathbf{\Psi}(1)^{-1}\right)\left(\mathbf{\Psi}(\pi)\mathbf{\Psi}(1)^{-1}\right)^{\prime}\bigg{)}^{-1}$
$\displaystyle\times\bigg{(}B(\pi)-\mathbf{\Psi}(\pi)\mathbf{\Psi}(1)^{-1}B(1)\bigg{)}.$

where $B(.)$ is a standard Brownian motion.

Using the asymptotic results given by Lemma 1, we simplify expression
$\mathbf{\Psi}(\pi)\mathbf{\Psi}(1)^{-1}$ to

$\mathbf{R}(\pi):=\mathbf{\Psi}(\pi)\mathbf{\Psi}(1)^{-1}=\begin{cases}\displaystyle\left(\pi\mathbf{\Omega}_{xx}+\int_{0}^{\pi}\underline{\mathbf{B}}dB^{\prime}\right)\left(\mathbf{\Omega}_{xx}+\int_{0}^{1}\underline{\mathbf{B}}dB^{\prime}\right)^{-1}&,\text{if}\
\gamma_{x}>1\\\ \\\
\displaystyle\left(\pi\mathbf{\Omega}_{xx}+\int_{0}^{\pi}\underline{\mathbf{J}}_{C}dJ_{C}^{\prime}\right)\left(\mathbf{\Omega}_{xx}+\int_{0}^{1}\underline{\mathbf{J}}_{C}dJ_{C}^{\prime}\right)^{-1}&,\text{if}\
\gamma_{x}=1\\\ \\\
\displaystyle\pi\mathbf{I}_{p}&,\text{otherwise}.\end{cases}$ (A.26)

Therefore, by denoting the Brownian functional above with
$\mathbf{N}(\pi)=\bigg{(}B(\pi)-\mathbf{R}(\pi)B(1)\bigg{)}$ and
$\mathbf{M}(\pi)=\bigg{(}\pi\left(\mathbf{I}_{p}-\mathbf{R}(\pi)\right)\left(\mathbf{I}_{p}-\mathbf{R}(\pi)\right)^{\prime}+(1-\pi)\mathbf{R}(\pi)\mathbf{R}(\pi)^{\prime}\bigg{)}$,
then we obtain

$\displaystyle\widetilde{\mathcal{W}}^{IVX}_{\beta}(t)\Rightarrow\underset{\pi\in[\pi_{1},\pi_{2}]}{\text{sup}}\bigg{\\{}\mathbf{N}(\pi)^{\prime}\mathbf{M}(\pi)^{-1}\mathbf{N}(\pi)\bigg{\\}}$
(A.27)

∎

Proof of Corollary 1 When $\gamma_{x}<1$, we have that
$\mathbf{R}(\pi)=\pi\mathbf{I}_{p}$ and thus
$\mathbf{M}(\pi)=\pi\left(\mathbf{I}_{p}-\pi\mathbf{I}_{p}\right)\left(\mathbf{I}_{p}-\pi\mathbf{I}_{p}\right)^{\prime}+(1-\pi)\pi\mathbf{I}_{p}\left(\pi\mathbf{I}_{p}\right)^{\prime}=\left[\pi(1-\pi)^{2}+(1-\pi)\pi^{2}\right]\mathbf{I}_{p}=\pi(1-\pi)\mathbf{I}_{p}$.
Hence, in this case, the limiting distribution in Theorem 2 will reduce to

$\displaystyle\widetilde{\mathcal{W}}^{IVX}_{\beta}(t)\Rightarrow\underset{\pi\in[\pi_{1},\pi_{2}]}{\text{sup}}\frac{\big{[}B(\pi)-\pi
B(1)\big{]}^{\prime}\big{[}B(\pi)-\pi B(1)\big{]}}{\pi(1-\pi)}.$ (A.28)

Thus, the above proof demonstrates that when we have predictors generated via
the LUR specification with $\gamma_{x}<1$ and testing the joint null
hypothesis of no structural break and no predictability then the limiting
distribution of the IVX Wald statistic converges to a NBB similar to the
result of Andrews, (1993) in the case of linear regression models.

Proof of Corollary 2. (i). Using the above notations, it’s straightforward to
obtain that under the null hypothesis $\beta_{1}=\beta_{2}=0$, we obtain that

$\displaystyle\mathcal{W}_{T}^{IVX}=\widetilde{\beta}^{IVX^{\prime}}\widetilde{\mathbf{Q}}_{\mathcal{R}}^{-1}\widetilde{\beta}^{IVX}\Rightarrow
U(1)^{\prime}\left[\mathbf{\Psi}(1)^{-1^{\prime}}\mathbf{\Sigma}\mathbf{\Psi}(1)^{-1}\right]^{-1}U(1).$
(A.29)

Combining this result with Theorem 2, we can deduce that

$\displaystyle\mathcal{W}_{\beta}^{joint}$
$\displaystyle=\mathcal{W}_{T}^{IVX}+\mathcal{W}_{\beta}^{IVX}(t)$
$\displaystyle\Rightarrow\begin{pmatrix}\mathbf{\Psi}(1)^{-1}U(1)\\\ \\\
U(\tau)-\mathbf{\Psi}(\pi)\mathbf{\Psi}(1)^{-1}U(1)\end{pmatrix}^{\prime}\begin{pmatrix}\mathbf{\Psi}(1)^{-1^{\prime}}\mathbf{\Sigma}\mathbf{\Psi}(1)^{-1}\
&\ \mathbf{0}_{p\times p}\\\ \\\ \mathbf{0}_{p\times
p}&\mathbf{\Delta}_{p\times
p}\end{pmatrix}^{-1}\begin{pmatrix}\mathbf{\Psi}(1)^{-1}U(1)\\\ \\\
U(\pi)-\mathbf{\Psi}(\tau)\mathbf{\Psi}(1)^{-1}U(1)\end{pmatrix}$

with

$\displaystyle\mathbf{\Delta}:=\pi\big{(}\mathbf{I}_{p}-\mathbf{\Psi}(\pi)\mathbf{\Psi}(1)^{-1}\big{)}\mathbf{\Sigma}\big{(}\mathbf{I}_{p}-\mathbf{\Psi}(\pi)\mathbf{\Psi}(1)^{-1}\big{)}^{\prime}+(1-\pi)\left(\mathbf{I}_{p}-\mathbf{\Psi}(\pi)\mathbf{\Psi}(1)^{-1}\right)\mathbf{\Sigma}\left(\mathbf{I}_{p}-\mathbf{\Psi}(\pi)\mathbf{\Psi}(1)^{-1}\right)^{\prime}$

Therefore, we obtain that

$\displaystyle\mathcal{W}_{\beta}^{joint}$
$\displaystyle=\mathcal{W}_{T}^{IVX}+\mathcal{W}_{\beta}^{IVX}(t)$
$\displaystyle\Rightarrow\begin{pmatrix}\mathbf{\Psi}(1)^{-1}B(1)\\\ \\\
B(\pi)-\mathbf{\Psi}(\pi)\mathbf{\Psi}(1)^{-1}B(1)\end{pmatrix}^{\prime}\begin{pmatrix}\mathbf{\Psi}(1)^{-1^{\prime}}\mathbf{\Psi}(1)^{-1}\
&\ \mathbf{0}_{p\times p}\\\ \\\ \mathbf{0}_{p\times
p}&\widetilde{\mathbf{\Delta}}_{p\times
p}\end{pmatrix}^{-1}\begin{pmatrix}\mathbf{\Psi}(1)^{-1}B(1)\\\ \\\
B(\pi)-\mathbf{\Psi}(\pi)\mathbf{\Psi}(1)^{-1}B(1)\end{pmatrix}$

with

$\displaystyle\widetilde{\mathbf{\Delta}}:=\pi\big{(}\mathbf{I}_{p}-\mathbf{\Psi}(\pi)\mathbf{\Psi}(1)^{-1}\big{)}\big{(}\mathbf{I}_{p}-\mathbf{\Psi}(\pi)\mathbf{\Psi}(1)^{-1}\big{)}^{\prime}+(1-\pi)\left(\mathbf{I}_{p}-\mathbf{\Psi}(\pi)\mathbf{\Psi}(1)^{-1}\right)\left(\mathbf{I}_{p}-\mathbf{\Psi}(\pi)\mathbf{\Psi}(1)^{-1}\right)^{\prime}$

which shows that

$\displaystyle\mathcal{W}_{\beta}^{joint}$
$\displaystyle=\mathcal{W}_{T}^{IVX}+\mathcal{W}_{\beta}^{IVX}(t)$
$\displaystyle\Rightarrow\begin{pmatrix}B(1)\\\ \\\
B(\pi)-\mathbf{R}(\pi)B(1)\end{pmatrix}^{\prime}\begin{pmatrix}\mathbf{I}_{p}\
&\ \mathbf{0}_{p\times p}\\\ \\\ \mathbf{0}_{p\times
p}&\mathbf{M}(\pi)\end{pmatrix}^{-1}\begin{pmatrix}B(1)\\\ \\\
B(\pi)-\mathbf{R}(\pi)B(1)\end{pmatrix}$ $\displaystyle\equiv
B(1)B(1)^{\prime}+\mathbf{N}(\pi)^{\prime}\mathbf{M}(\pi)^{-1}\mathbf{N}(\pi)$

Since the first component of the above decomposition is independent of $\pi$,
then by the Continuous Mapping Theorem, we conclude that

$\displaystyle\widetilde{\mathcal{W}}_{\beta}^{joint}=\underset{\pi\in[\pi_{1},\pi_{2}]}{\text{sup}}\bigg{\\{}\mathcal{W}_{T}^{IVX}+\mathcal{W}_{\beta}^{IVX}(t)\bigg{\\}}\Rightarrow
B(1)B(1)^{\prime}+\underset{\pi\in[\pi_{1},\pi_{2}]}{\text{sup}}\bigg{\\{}\mathbf{N}(\pi)^{\prime}\mathbf{M}(\pi)^{-1}\mathbf{N}(\pi)\bigg{\\}}.$

(ii). In case that $\gamma_{x}<1$, by Corollary 1 it holds that the second
component of the $\mathcal{W}_{\beta}$ test statistic reduces to a functional
of a standard Brownian bridge $B(\pi)-\pi B(1)$. Then, since both $B(1)$ and
$B(\pi)-\pi B(1)$ are Gaussian processes which implies that
Cov$\big{(}B(1),B(\pi)-\pi B(1)\big{)}=\pi-\pi=0$. Therefore, these two
stochastic quantities are independent of each other313131To validate the
asymptotic independence of the two BM functionals, we can apply properties of
the BM and prove that the covariance of the two terms is zero, which ensures
independence.. Hence, the proof of the statement follows.

Proof of Proposition 1. Under the null hypothesis,
$\mathbb{H}_{0}:\alpha_{1}=\alpha_{2}$ and $\beta_{1}=\beta_{2}=\beta$, there
is no break in the model intercept. In this section we consider in more
details the estimator of the model intercept based on the IVX instrumentation
before proving the asymptotic distribution given by Proposition 1. In
particular, we propose to estimate the model intercept using the generated
instrument instead of the predictor, we refer to this estimate as
$\alpha^{IVZ}$ and the econometric intuition is explained below.

Based on the IVX estimation procedure, the corresponding IVX estimate for the
model intercept is given by
$\hat{\alpha}=\bar{y}_{T}-\widetilde{\beta}^{IVX}\bar{x}_{T-1}$. However, we
notice that due to the presence of the predictors, then the limit theory for
the estimate of the model intercept will vary with the degree of persistence
of predictors. We can see this below

$\displaystyle\sqrt{T}\left(\hat{\alpha}^{IVX}-\alpha\right)=\sqrt{T}\bar{u}_{T}-\left[T^{\frac{1+\gamma_{x}\wedge\delta_{z}}{2}}\left(\widetilde{\beta}^{IVX}-\beta\right)\right]\left[T^{-\frac{\gamma_{x}\wedge\delta_{z}}{2}}\bar{x}_{T-1}\right],$
(A.30)

where $\bar{u}_{T}=\frac{1}{T}\sum_{t=1}^{T}u_{t}$. Notice that both
$\sqrt{T}\bar{u}_{T}$ and
$T^{\frac{1+\gamma_{x}\wedge\delta_{z}}{2}}\left(\widetilde{\beta}-\beta\right)$
are both $\mathcal{O}_{p}(1)$, while the order of convergence of the last term
depends on the persistence level of predictors. We have the following
convergence rates

$\displaystyle\sum_{t=1}^{T}x_{t-1}=\begin{cases}\mathcal{O}_{p}\left(T^{-1/2}\right)&,\text{if}\
\gamma_{x}=1\\\ \\\ \mathcal{O}_{p}\left(T^{1/2+\gamma_{x}}\right)&,\text{if}\
0<\gamma_{x}<1\\\ \\\
\mathcal{O}_{p}\left(T^{3/2}\right)&,\gamma_{x}>1.\end{cases}$ (A.31)

We can observe that the term
$T^{-\frac{\gamma_{x}\wedge\delta_{z}}{2}}\bar{x}_{T-1}$ will dominate in the
limit in the case that $\gamma_{x}>\left(\delta_{z}+1\right)/2$ and vanish if
the reverse holds. In the case that $\gamma_{x}=\left(\delta_{z}+1\right)/2$,
all three terms appear in the asymptotic distribution which will depend on an
unknown parameter $\beta$. To simplify the asymptotic theory we need to derive
we estimate the model intercept based on the generated endogenous instrument,
that is, $\alpha^{IVZ}=\bar{y}_{T}-\widetilde{\beta}^{IVX}\bar{z}_{T-1}$,
where $\bar{z}_{T-1}=\frac{1}{T-1}\sum_{j=2}^{T}\tilde{z}_{j-1}$. The
advantage of the IVZ model estimate is that the persistence level of the
instrument $\widetilde{z}_{t}$ is controlled by the choice of the tuning
parameters $\delta_{z}$ and $c_{z}$ therefore the abstract degree of
persistence is filtered out.

Moreover, the particular choice of the estimate for the model intercept works
as a power enhancement mechanism against nonzero $\beta$ values, as seen below

$\displaystyle\left(\hat{\alpha}^{IVZ}-\alpha\right)=\bar{u}_{T}-\left(\widetilde{\beta}^{IVX}-\beta\right)\bar{z}_{T-1}+\beta\left(\bar{x}_{T-1}-\bar{z}_{T-1}\right),$
(A.32)

The order of the first term is $\mathcal{O}_{p}\left(T^{-1/2}\right)$.
Moreover, we can show that the second term is asymptotically dominated by the
first term, while the third term does not converge unless the parameter space
for $\beta$ is within the neighbourhood of zero. Therefore, power against the
presence of predictability when a model intercept is included in the model is
achieved when we control the convergence rate of the third term above.

Next, we show that the second term of (A.32) is asymptotically dominated by
the first term by expanding further the expression

$\displaystyle\sqrt{T}\left(\widetilde{\beta}^{\text{IVX}}-\beta\right)\bar{z}_{T-1}\equiv\left[T^{\frac{1+\gamma_{x}\wedge\delta_{z}}{2}}\left(\widetilde{\beta}^{\text{IVX}}-\beta\right)\right].\left[T^{-\left(1+\frac{\gamma_{x}\wedge\delta_{z}}{2}\right)}\sum_{t=1}^{T}\widetilde{z}_{t-1}\right].$
(A.33)

The first term above is $\mathcal{O}_{p}(1)$ due to the convergence property
of the IVX estimator to the mixed Gaussian distribution as established by
Phillips and Magdalinos, (2009). To establish the order of the second term we
consider the convergence rate of the term $\sum_{t=1}^{T}\widetilde{z}_{t-1}$.
For instance, by Lemma A2 in the Online Appendix of KMS we have that

$\displaystyle\sum_{t=1}^{T}\widetilde{z}_{t-1}=\begin{cases}\mathcal{O}_{p}\left(T^{\frac{\gamma_{x}\wedge
1}{2}+\delta_{z}}\right)&,\text{if}\ \delta_{z}<\gamma_{z}\\\ \\\
\mathcal{O}_{p}\left(T^{\gamma_{x}+\frac{\delta_{z}}{2}}\right)&,\text{if}\
0<\gamma_{x}<\delta_{z}.\end{cases}$ (A.34)

Hence, we have that the second term of (A.32) will be
$\mathcal{O}_{p}\left(T^{\frac{\gamma_{x}\wedge 1+\delta_{z}}{2}+1}\right)$.
Since $\gamma_{x}\wedge 1+\delta_{z}<2$, this implies that the order of the
second term is $o_{p}(1)$. Therefore, we prove that the second term of
$\left(\hat{\alpha}^{IVZ}-\alpha\right)$ is asymptotically dominated by the
first term.

For the third term of (A.32) we aim to show that is asymptotically dominated
by the first term when $\beta$ is small enough. Thus, we rewrite with
$\sqrt{T}\beta\left(\bar{x}_{T-1}-\bar{z}_{T-1}\right)=\beta\frac{1}{\sqrt{T}}\sum_{t=1}^{T}\left(x_{t-1}-\widetilde{z}_{t-1}\right)$.
Using the representation formula of the instrument given by expression (23) of
Phillips and Magdalinos, (2009), we find an equivalent expression, that is,
$\beta\frac{1}{\sqrt{T}}\sum_{t=1}^{T}\left(x_{t-1}-\widetilde{z}_{t-1}\right)=-\beta
C_{z}T^{\left(\frac{1}{2}+\delta_{z}\right)}\sum_{t=1}^{T}\psi_{t-1}$, where
$\psi_{t-1}=\sum_{j=1}^{t}R_{T}^{t-j}x_{j-1}$. Hence, if we set for simplicity
$\mathbf{C}_{z}$ = $\mathbf{I}_{p}$ indicating a common degree of persistence
across the instruments, we then obtain the following probability bound for
this expression

$\displaystyle\left\lVert\sqrt{T}\beta\left(\bar{x}_{T-1}-\bar{z}_{T-1}\right)\right\rVert$
$\displaystyle\leq\left\lVert\beta\right\rVert
T^{-\left(\frac{1}{2}+\delta_{z}\right)}\sum_{t=1}^{T}\left\lVert\psi_{t-1}\right\rVert$
$\displaystyle\leq\left\lVert\beta\right\rVert
T^{-\left(\frac{1}{2}+\delta_{z}\right)}T\underset{2\leq t\leq
T}{\text{sup}}\left\lVert\psi_{t-1}\right\rVert$

$\displaystyle=\left\lVert\beta\right\rVert
T^{\left(\frac{1}{2}-\delta_{z}\right)}\underset{2\leq t\leq
T}{\text{sup}}\left\lVert\psi_{t-1}\right\rVert$
$\displaystyle\leq\left\lVert\beta\right\rVert
T^{\left(\frac{1}{2}-\delta_{z}\right)}\mathcal{O}_{p}\left(T^{\frac{\gamma_{x}}{2}+\delta_{z}}\right)$
$\displaystyle=\left\lVert\beta\right\rVert\mathcal{O}_{p}\left(T^{\frac{\gamma_{x}}{2}+\delta_{z}}\right),$

Notice that the above result is justified due to the uniform bound of
$\left\lVert\psi_{t-1}\right\rVert$, which is shown to be
$\mathcal{O}_{p}\left(T^{\frac{\gamma_{x}}{2}+\delta_{z}}\right)$ by PM. Thus,
if $\beta=o_{p}\left(T^{-\frac{1+\gamma_{x}}{2}}\right)$, then
$\left\lVert\sqrt{T}\beta\left(\bar{x}_{T-1}-\bar{z}_{T-1}\right)\right\rVert=o_{p}(1)$
and the third component of $\left(\tilde{\alpha}^{IVZ}-\alpha\right)$ will
also be asymptotically dominated. Therefore, when $\beta$ is a non-zero
constant, this term will dominate, and this is the reason we obtain non-zero
local power under the alternative hypothesis of predictability.

Proof of Proposition 2.

For the proof of Proposition 2, we have already proved most of the required
results.The only step we need to additionally show is that the test statistics
$\mathcal{W}_{a}(t)$ and $\mathcal{W}_{b}(t)$ are asymptotically independent
of each other. Notice that since the test statistic $\mathcal{W}_{a}(t)$ is
driven by $\sum_{t=1}^{\lfloor Tr\rfloor}u_{t}$, while the test statistic
$\mathcal{W}_{b}(t)$ is driven by $\sum_{t=1}^{\lfloor
Tr\rfloor}\widetilde{z}_{t-1}u_{t}$. These two partial sums (invariance
principles) have a joint weakly convergence to two independent Brownian
motions, as shown by Proposition A1 in Phillips and Magdalinos, (2009). Hence,
the asymptotic independence is guaranteed. Then, convergence of the test
statistic $\mathcal{W}_{\alpha\beta}(t)$ follows by an application of the
continuous mapping theorem.

##  Appendix B Critical Values and Implementation

###  B.1 Critical Values

The correct use of critical values close to the true asymptotic distribution
of the test statistic under examination is crucial in Monte Carlo simulations.
This allows to correctly identify the existence of size distortions. Notice
that for deciding whether to accept or reject the null hypothesis, one can use
already tabulated critical values, only in those cases that the asymptotic
distribution of the test statistic under consideration has an approximate form
to Brownian functional (such as NBB). Moreover, we apply bootstrap
methodologies (such as the fixed regressor bootstrap) to ensure that we obtain
robust performance of the empirical size and power of the proposed tests.

###  B.2 Software Implementation

The Monte Carlo simulation study of the paper was implemented via the
Statistical Package R (Version 3.5.1). More specifically, to validate the code
for testing for both predictability and structural break, using the supremum
Wald type statistics we utilize both the Matlab code of KMS as well as the IVX
library written by Vasilopoulos and Pavlidis, (2019) which provide
implementations of the IVX instrumentation procedure. Moreover, to reduce the
execution time we utilize parallel programming techniques as well as related R
packages for this purpose, such as the Rcpp library and the Rccp Armadillo
linear algebra library. Then, the critical values and empirical size and power
of the tests were obtained by simulating the asymptotic distributions based on
10,000 Monte Carlo replications using the Iridis4 High Performance Computing
Facility of the University of Southampton.

###  B.3 Bootstrap procedure

For robustness of the procedure we use a bootstrap algorithm in order to
obtain corresponding critical values. The Algorithm for the bootstrapped
critical values that we follow in this paper is proposed by Xu, (2020), even
though we do not consider the long-horizon predictability component of the
particular procedure. We also consider the fixed regressor bootstrap as
described by Georgiev et al., (2018) and also in Hansen, 2000b .

##### Algorithm 1.

red(to modify) Two-sided $100\alpha\%$ Bootstrap Implied Test for
$\mathbb{H}_{0}:R\theta=r_{0}$

Consider the predictive regression with multiple predictors given by

$\displaystyle y_{t}$ $\displaystyle=\mu_{y}+x_{t-1}^{\prime}\beta+u_{yt},$
(B.1) $\displaystyle x_{t}$ $\displaystyle=\mu_{x}+Rx_{t-1}+u_{xt},$ (B.2)

where $R$ is a $k\times k$ dimensional coefficient matrix and the error
sequence $u_{t}=\left(u_{yt},u_{xt}^{\prime}\right)^{\prime}$ is $(k\times
1)-$dimensional martingale difference sequence with variance-covariance matrix

$\displaystyle\text{Var}(u_{t})=\Sigma=\begin{bmatrix}\sigma^{2}_{y}&\
\sigma_{xy}^{\prime}\\\ \sigma_{xy}&\ \Sigma_{xx}^{\prime}\end{bmatrix}$ (B.3)

  * Step 1.

Estimation Step

    * 1.1.

Run OLS regressions of (B.1) and (B.2) to obtain the OLS estimates of the
vector of coefficients
$\left\\{\hat{\mu}_{y},\hat{\beta},\hat{\mu}_{x},\hat{R}\right\\}$ and the
sequence of residuals
$\hat{u}_{t}=\left(\hat{u}_{yt},\hat{u}_{xt}^{\prime}\right)^{\prime}$.

    * 1.2.

Compute the Wald statistic $\mathcal{W}_{IM}$.

  * Step 2.

Generate the bootstrap sample
$\left\\{y_{t}^{*},x_{t}^{*},t=1,...,n\right\\}$.

    * 2.1.

Generate $u_{t}^{*}=\hat{e}_{t}e_{t}$, where $e_{t}$ is a random number
generated from the standard normal distribution. Write
$u_{t}^{*}=\left(u_{yt}^{*},u_{xt}^{*\prime}\right)^{\prime}$.

    * 2.2.

Generate $x_{t}^{*}$ such that
$x_{t}^{*}=\hat{\mu}_{x}+\hat{R}_{bc}x_{t-1}^{*}+u_{xt}^{*}$, where
$\hat{R}_{bc}$ is the bias-corrected estimator of $R$, defined as below

$\displaystyle\hat{R}_{bc}=\hat{R}+\hat{\Sigma}_{xx}\bigg{[}\big{(}I_{k}-\widehat{R}^{\prime}\big{)}^{-1}+\widehat{R}^{\prime}\big{(}I_{k}-\widehat{R}^{\prime
2}\big{)}^{-1}+\sum_{j=1}^{k}\widehat{\lambda}_{j}\big{(}I_{k}-\widehat{\lambda}_{j}\widehat{R}^{\prime}\big{)}^{-1}\bigg{]}\left(\sum_{t=2}^{n}\tilde{x}_{t-1}\tilde{x}_{t-1}^{\prime}\right),$

where $\widehat{\lambda}_{j}^{\prime}$s for $j=1,...,k$ are $k$ eigenvalues of
$\widehat{R}^{\prime}$. We set $x_{1}^{*}=x_{1}$.

    * 2.3.

Generate $y_{t}^{*}$ such that
$y_{t}^{*}=\hat{\mu}_{y}+x_{t-1}^{*\prime}\hat{\beta}+u_{yt}^{*}$. We set
$y_{1}^{*}=y_{1}$.

    * 2.4.

Repeat Steps 2.1 to 2.3 $B$ times (e.g., $B=1000$).

  * Step 3.

For each bootstrap sample, calculate the bootstrap Wald statistic

$\displaystyle\mathcal{W}^{*}_{IM}=\left(\mathcal{R}\hat{\theta}^{*}_{h}-\mathcal{R}\hat{\theta}_{h}\right)^{\prime}\left[\mathcal{R}\widehat{q}^{*}\left(\widehat{\Sigma}^{*}\otimes\left(\sum_{t=2}^{n}\widetilde{x}^{*}_{t-1}\widetilde{x}^{*\prime}_{t-1}\right)^{-1}\right)\widehat{q}^{*\prime}\mathcal{R}^{\prime}\right]\left(\mathcal{R}\hat{\theta}^{*}_{h}-\mathcal{R}\hat{\theta}_{h}\right),$

where $\hat{\theta}^{*}_{h}$, $\widehat{q}^{*}$, $\widehat{\Sigma}^{*}$ are
computed like $\hat{\theta}_{h}$, $\widehat{q}$, $\widehat{\Sigma}$ except
that the bootstrap sample $\left\\{y_{t}^{*},x_{t}^{*}\right\\}$ is replaced
by the original sample $\left\\{y_{t},x_{t}\right\\}$.

  * Step 4.

Use the $100(1-\alpha)\%$ quantile $\tau_{\mathcal{W}_{IM}^{*}}(1-\alpha)$ of
$\mathcal{W}_{IM}^{*}$ over $B$ bootstrap replications as the critical value,
that is, the hypothesis $H_{0}$ is rejected at the significance level $\alpha$
if $\mathcal{W}_{IM}>\tau_{\mathcal{W}_{IM}^{*}}(1-\alpha)$.

###  B.4 Structural Break Tests

For the empirical application of the paper, we consider a set of additional
structural break tests as robustness checks. Moreover, the particular
structural break tests can be employed to identify the location of the break-
point.

##  Appendix C Tables and Figures

###  C.1 Empirical Results

Table C1: Predictability Tests for the Equity Premium of univariate predictive regressions Predictor | $\widehat{\beta}^{\text{OLS}}$ | $t^{\text{HAC}}$ | $R^{2}$ | $\widehat{\beta}^{\text{IVX}}$ | $\mathcal{W}^{\text{IVX}}$ | $sup\mathcal{W}^{\text{OLS}}$ | $sup\mathcal{W}^{\text{IVX}}$
---|---|---|---|---|---|---|---
Panel A: | 1946Q1 - 2019Q4
Dividend payout ratio | 0.0056 | 0.3898 | 0.0003 |  |  |  |
Long-term yield | -0.0824 | -1.5958 | 0.0032 |  |  |  |
Dividend yield | 0.0150 | 1.9880 | 0.0047 |  |  |  |
Dividend-price-ratio | 0.0142 | 1.8741 | 0.0042 |  |  |  |
T-bill rate | -0.1043 | -2.1461∗ | 0.0060 |  |  |  |
Earnings-price-ratio | 0.0115 | 1.2885 | 0.0028 |  |  |  |
Book-to-market ratio | 0.0042 | 0.6590 | 0.0006 |  |  |  |
Default yield spread | 0.0150 | 1.9880∗ | 0.0047 |  |  |  |
Net equity expansion | -0.0501 | -0.6149 | 0.0005 |  |  |  |
Term spread | 0.1813 | 1.6834 | 0.0035 |  |  |  |
Inflation rate | -0.9075 | -2.6520∗∗ | 0.0096 |  |  |  |
Panel B: | 1990Q1 - 2019Q4
Dividend payout ratio | 0.0028 | 0.1597 | 0.0001 |  |  |  |
Long-term yield | -0.0883 | -0.8409 | 0.0017 |  |  |  |
Dividend yield | 0.0354 | 1.6497 | 0.0102 |  |  |  |
Dividend-price-ratio | 0.0351 | 1.6389 | 0.0101 |  |  |  |
T-bill rate | -0.0393 | -0.3885 | 0.0005 |  |  |  |
Earnings-price-ratio | 0.0173 | 0.8403 | 0.0041 |  |  |  |
Book-to-market ratio | 0.0317 | 1.0439 | 0.0041 |  |  |  |
Default yield spread | 0.5275 | 0.5506 | 0.0025 |  |  |  |
Net equity expansion | 0.1188 | 0.9655 | 0.0037 |  |  |  |
Term spread | -0.0727 | -0.4442 | 0.0005 |  |  |  |
Inflation rate | 0.2005 | 0.2831 | 0.0003 |  |  |  |

Table C1 presents simple predictability tests of the null hypothesis $\beta=0$
and model estimates based on the predictive regression with a single predictor
given by $y_{t}=\alpha+\beta x_{t-1}+u_{t}$, applied to all individual
predictors. For the t-tests we consider a t-ratio based on the HAC (Newey-
West) covariance estimator. For all test statistics, we denote the rejection
probabilities under the null hypothesis of no predictability at significance
levels $1\%(^{***})$, $5\%(^{**})$, and $10\%(^{*})$, respectively.

Table C2: Structural Break Tests for predictors Panel A: | 1990Q1 - 2019Q4
---|---
Predictor |  | $\widehat{\rho}$ | $R^{2}$ | BP Seq F-test | Max LR F-test | Exp LR F-test
Dividend payout ratio |  | 0.7477 | 0.5592 | 4.0061 | 49.4450∗∗∗ | 19.1970∗∗∗
Long-term yield |  | 0.0429 | 0.0046 | 3.4002 | 2.3062 | 0.3676
Dividend yield |  | 0.0424 | 0.0001 | 3.5761 | 7.2497 | 1.4443
Dividend-price-ratio |  | 0.0511 | 0.0008 | 3.6305 | 7.5262 | 1.5400
T-bill rate |  | 0.4734 | 0.2169 | 2.6279 | 5.4363 | 0.6033
Earnings-price-ratio |  | 0.4385 | 0.1922 | 1.6033 | 14.6679∗∗∗ | 4.6243∗∗∗
Book-to-market ratio |  | -0.0439 | 0.0012 | 5.5639 | 3.1215 | 0.1470
Default yield spread |  | 0.4616 | 0.2131 | 3.1040 | 7.4651 | 1.9635∗∗
Net equity expansion |  | 0.2440 | 0.0596 | 0.9001 | 2.5452 | 0.2197
Term spread |  | 0.1035 | 0.0107 | 1.6129 | 2.1243 | 0.2136
Inflation rate |  | -0.0979 | 0.0096 | 3.2371 | 5.1032 | 0.9540
Predictor | $\widehat{\mu}$ | $\widehat{\rho}$ | $R^{2}$ | BP Seq F-test | Max LR F-test | Exp LR F-test
Dividend payout ratio | -0.0001 | 0.7476 | 0.5592 | 26.5211∗∗∗ | 26.5211∗∗∗ | 7.7360∗∗∗
Long-term yield | -0.0002 | 0.0368 | 0.0014 | 1.7122 | 1.1633 | 0.1840
Dividend yield | -0.0007 | 0.0408 | 0.0017 | 3.4744 | 5.0447 | 1.1532
Dividend-price-ratio | -0.0007 | 0.0495 | 0.0024 | 3.3367 | 5.0907 | 1.1876
T-bill rate | -0.0001 | 0.4684 | 0.2197 | 1.3134 | 2.6698 | 0.4553
Earnings-price-ratio | -0.0003 | 0.4384 | 0.1922 | 0.7989 | 7.5124∗∗∗ | 2.2289
Book-to-market ratio | -0.0006 | -0.0447 | 0.0020 | 2.7760 | 2.0446 | 0.2170
Default yield spread | 0.0000 | 0.4616 | 0.2131 | 1.6369 | 3.7191 | 0.7621
Net equity expansion | 0.0000 | 0.2440 | 0.0596 | 0.8209 | 1.6173 | 0.1941
Term spread | 0.0000 | 0.1034 | 0.0107 | 1.3025 | 1.7316 | 0.2593
Inflation rate | 0.0000 | -0.0979 | 0.0097 | 1.6209 | 2.5489 | 0.4311

Table C2 presents structural break tests and model estimates based on an AR(1)
model applied to all individual stationary predictors (based on first
differences, i.e., $\Delta x=x_{t}-x_{t-1}$) for Panel A with sampling period:
1990Q1 - 2019Q4. The covariance matrix for the AR(1) model is constructed
using the HAC (Newey-West) estimator. We consider the following structural
break tests323232Further details regarding the specification and testing
algorithms can be found in the User’s Guide of Eviews under the section titled
”Stability Diagnostics”, see https://www.eviews.com/help/ (i) Bai-Perron
F-test based on the sequential break detection algorithm of Bai and Perron,
(2003) with $\epsilon=0.15$; (ii) Maximum LR F-statistic with $\epsilon=0.15$;
and (iii) Exp LR F-statistic with $\epsilon=0.15$. The last two structural
break tests represent the Andrews unknown breakpoint tests (see, Andrews,
(1993)). Notice also for the BP statistic we assume a common data distribution
across the blocks to ensure consistent estimation of the variance. For all
test statistics, we denote the rejection probabilities under the null
hypothesis of no structural break at significance levels $1\%(^{***})$,
$5\%(^{**})$, and $10\%(^{*})$, respectively.


# Critical behavior of the stochastic SIR model on random bond-diluted
lattices

Carlos Handrey A. Ferraz<EMAIL_ADDRESS>José Luiz S. Lima  Exact and
Natural Sciences Center, Universidade Federal Rural do Semi-Árido-UFERSA, PO
Box 0137, CEP 59625-900, Mossoró, RN, Brazil

###### Abstract

In this paper, we investigate the impact of bond-dilution disorder on the
critical behavior of the stochastic SIR model. Monte Carlo simulations were
conducted using square lattices with first- and second-nearest neighbor
interactions. Quenched bond-diluted lattice disorder was introduced into the
systems, allowing them to evolve over time. By employing percolation theory
and finite-size scaling analysis, we estimate both the critical threshold and
leading critical exponent ratios of the model for different bond-dilution
rates ($p$). An examination of the average size of the percolating cluster and
the size distribution of non-percolating clusters of recovered individuals was
performed to ascertain the universality class of the model. The simulation
results strongly indicate that the present model belongs to a new universality
class distinct from that of 2D dynamical percolation, depending on the
specific $p$ value under consideration.

###### keywords:

Stochastic SIR model , bond-diluted lattices , percolating cluster , critical
exponent ratios , Monte Carlo simulation

##  1 Introduction

In epidemiology, understanding the spread of infectious diseases is crucial
for disease control and public health management. Bond-diluted lattices offer
a platform to study the impact of network structure on disease transmission
dynamics. By examining how diseases propagate on networks with partial
connections, it is possible to develop more effective strategies for disease
prevention, intervention, and control. A fundamental system of ordinary
differential equations was initially developed for modeling the dynamics of
infectious diseases within populations [1]. This system categorized
individuals into three basic classes: susceptible (S), infected (I), and
recovered (R). This gave rise to the SIR model, named after its component
classes. Over time, this model has undergone numerous adaptations and
extensions, such as the SIS, SIRS and SEIR models, each tailored to capture
different nuances of disease transmission within populations.

The simplicity of SIR model lies in its assumption that infected individuals
either gain permanent immunity or succumb to the disease, thereby exiting the
transmission dynamics. This characteristic renders the model particularly
suitable for simulating outbreaks of influenza[2], measles [3, 4], varicella
[5], AIDS [6, 7], COVID-19 [8, 9, 10], and many others. This model, rooted in
both collective dynamics [11, 12] and complex systems [13], serve as
invaluable tool for assessing disease propagation rates and evaluating the
efficacy of public health interventions. It has been shown that the SIR model
shares a universality class with dynamic percolation [14, 15, 16, 17, 18].

Lattice models excel in capturing various aspects of infection spread across
spatially extended landscapes, where wave-like evolution plays a key role.
This work focuses on a stochastic lattice-gas variant of the SIR model with
asynchronous site updates, a framework previously employed in analyzing
various population dynamics [19, 20]. In both synchronous and asynchronous
implementations, the model exhibits a phase transition as control parameters
vary. This transition, typically second-order, delineates two distinct
regimes: an endemic phase where the population largely remains susceptible,
and an epidemic phase characterized by widespread disease transmission,
leading to a significant proportion of the population becoming infected,
recovering, or succumbing to the disease. At the critical transition point,
the system reaches an epidemic outbreak threshold. It is worth mentioning that
both stochastic and deterministic descriptions of the model can be connected
through mean-field approximation [21, 22], resulting in Langevin equations
associated with the Fokker–Planck equation.

In the present study, we will conduct an analysis focusing on two key aspects:
the average size of the percolating cluster and the size distribution of non-
percolating clusters among recovered individuals. To achieve this, we will
employ the Newman–Ziff algorithm [23], which facilitates the determination of
critical exponent ratios within the model.

The SIR stochastic lattice model behaves akin to an absorbing system,
characterized by an active phase featuring an infinite array of absorbing
configurations wherein the final state solely comprises recovered individuals.
However, our approach differs: we halt simulations upon confirming the
existence of a percolating cluster (spanning cluster) within the system,
effectively establishing a non-absorbing state for the model. This strategy
expedites data analysis without compromising the assurance of reaching the
system’s asymptotic limit.

Notably, in the critical regime of the system, only a single percolation
cluster is generated. Through a meticulous analysis of these clusters, we aim
to ascertain the universality class of the model. Furthermore, the network
topology serves as a crucial determinant directly influencing the dynamics of
involved processes [24, 25, 26]. Therefore, we will investigate how quenched
bond-diluted lattice disorder influence both dynamics and universality class
of this model.

In this context, bond-diluted lattices emerge from regular square lattices
($N=L\times L$) initially featuring both first- and second-nearest neighbor
interactions. In these lattices, each site is subject to a probability $p$ of
experiencing disconnection in its connections (dilution), with the requirement
that at least one connection remains intact, thereby preventing isolated
sites. These networks aim to emulate a more realistic population, reflecting a
diverse array of interconnected relationships among individuals. The present
work investigates several cases of varying $p$.

This study mainly aims to understand how the dilution factor affects the
critical behavior of the stochastic SIR model. Since the dilution probability
$p$ can be understood as a quenched topological disorder, we want to determine
if this kind of disorder is relevant to changing the universality class of the
model. Several criteria have been devised to anticipate whether quenched
disorder might alter the critical exponents of a given model, including the
well-known Harris criterion [27] and its refinement, the
Harris–Barghathi–Vojta criterion (HBV) [28]. The Harris criterion posits that
a second-order transition in a $d$-dimensional system, with an original
correlation length exponent $\nu$, remains stable against quenched spatial
disorder if $\nu>2/d$. Conversely, the HBV criterion suggests that quenched
topological disorder is inconsequential concerning phase transition stability
if the system satisfies $\nu>1/k$, where $k$ denotes the disorder decay
exponent measuring how swiftly coordination number fluctuations decrease with
increasing system length scale. However, it is worth noting that these
criteria have been observed to falter in certain cases [29].

The manuscript is organized as follows: In Section 2, we provide an overview
of the SIR stochastic lattice model. Section 3 details the Monte Carlo (MC)
simulation background and finite-size analysis. In Section 4, we present and
discuss the results. Finally, Section 5 we present our conclusions.

##  2 SIR stochastic lattice model

The SIR stochastic lattice model is formulated on a lattice consisting of $N$
sites, where each site can host a single individual characterized as
susceptible ($S$), infected ($I$), or recovered (immunized/deceased) ($R$).
This model encompasses two distinct processes: an auto-catalytic mechanism,
denoted by $I+S\rightarrow I+I$, and a spontaneous transition from infection
to recovery, represented by $I\rightarrow R$. At each time step, a site is
randomly chosen, and a specific set of dynamic rules is applied in the
following way:

  1. (i)

If the site is in the state $S$ and there is at least one neighboring site in
the sate $I$ then the site becomes $I$ with probability proportional to a
parameter $\mu$ and the number $z$ of neighboring sites, i.e., $\mu{\kern
1.0pt}m/z$, where $m$ is the number of neighboring $I$ sites.

  2. (ii)

If the site is in state $I$ it becomes $R$ spontaneously with probability
$\lambda$.

  3. (iii)

If the site is $R$ it remains unchanged.

At each site $i$ of a lattice we assign a stochastic variable $\sigma_{i}$
that can take the values 0, 1 or 2, according to whether the site is in the
state $S$, $I$, or $R$, respectively. Since the transitions between states in
this model are non-equilibrium ones, the allowed transitions of the state $i$
of a site are cyclic, this is, $0\rightarrow 1\rightarrow 2$. The
corresponding transition rate is represented by $w_{i}(\sigma)$ and describes
the transition $\sigma\rightarrow\sigma^{\prime}$ in which the whole
microscopic configuration (microstate)
$\sigma^{\prime}\equiv(\sigma_{1},\ldots,\sigma^{\prime}_{i},\ldots,\sigma_{N})$
differs from $\sigma$ only by the state of the $i$-th site. It is given by

$w_{i}(\sigma)=\frac{\mu}{z}\delta(\sigma_{i},0)\sum\limits_{j}{\delta(\sigma_{j},1)}+\lambda{\kern
1.0pt}\delta(\sigma_{i},1),$ (1)

where the summation runs over the nearest neighbors of site $i$ and
$\delta(x,y)$ denotes the Kronecker delta. The parameters $\mu$ and $\lambda$
are related to the subprocesses above described, and are chosen such that
$\mu+\lambda=1$.

The system evolves in time according to a master equation for the probability
distribution $W(\sigma,t)$ described by

$\frac{d}{{dt}}W(\sigma,t)=\sum\limits_{i}{\\{w_{i}(\bar{\sigma})}W(\bar{\sigma},t)-w_{i}(\sigma)W(\sigma,t)\\},$
(2)

where the microstate $\bar{\sigma}$ is obtained from $\sigma$ by an anticyclic
permutation ($2\rightarrow 1\rightarrow 0$) of the state of the site $i$.

##  3  Numerical simulation and finite-size scaling analysis

We can establish an asynchronous, non-absorbing SIR model by applying the
following kinetic Monte Carlo rules:

  1. (i)

First, we start with a single central infected site (seed) and the remaining
ones being all susceptible on a two-dimensional lattice in which each
individual of the population is attached to its respective lattice site. In
order to speed up the simulation we create two lists that are updated at each
algorithm step: a list of infected individuals (infected list) and a list of
recovered individuals (recovered list), which begins empty.

  2. (ii)

Next, we update the system state by randomly choosing an available infected
site from the infected list and proceed as follows:

    1. (a)

Generate a random number ($rn$) in the interval $(0,1)$. If $rn\leq\lambda$,
the infected site is removed from the infected list and placed in the
recovered list;

    2. (b)

Otherwise (if $rn>\lambda$), pick randomly one nearest neighbor of the
infected site and make it also infected provided that it is susceptible,
adding it to infected list.

  3. (iii)

Repeat asynchronously the step (ii) several times until either there is no
infected sites (endemic phase) or there is a percolating cluster of recovered
sites (non-absorbing epidemic phase).

Remarkably, research has shown that the SIR model on square lattices shares
the same universality class as dynamic percolation (DP). This enables us to
explore the phase transition occurring in the current non-absorbing SIR model
using percolation theory. Consequently, we can define the epidemic phase of
the model when a percolating cluster of recovery sites forms in the system and
the endemic phase when it does not. Such a graphical analogy has also been
applied to other compartmental models.

Following classical percolation theory, it is crucial to first determine the
cluster distribution of recovery sites, i.e., the number of clusters with $s$
recovery sites denoted as $n_{p}(s)$. This can be achieved using the
Newman–Ziff algorithm, which also features a built-in capability to identify
whether a percolating cluster has formed or not based on the considered
$\lambda$ value. It’s noteworthy that for systems with non-periodic
boundaries, like the ones under consideration here, the percolating cluster is
essentially a spanning cluster [23, 30].

Table 1:  Estimates of the epidemic threshold $\lambda_{c}$ and critical exponents ratios for each $p$ cases. Lattice | Epidemic threshold | $1/\nu$ | $\beta/\nu$ | $\gamma/\nu$
---|---|---|---|---
$p=0$ | $\lambda_{c}=0.275(0)$ | $0.716\pm 0.014$ | $0.112\pm 0.003$ | $1.776\pm 0.009$
$p=0.10$ | $\lambda_{c}=0.266(3)$ | $0.667\pm 0.022$ | $0.048\pm 0.006$ | $1.722\pm 0.009$
$p=0.20$ | $\lambda_{c}=0.258(1)$ | $0.691\pm 0.019$ | $0.054\pm 0.005$ | $1.744\pm 0.012$
$p=0.30$ | $\lambda_{c}=0.250(5)$ | $0.788\pm 0.011$ | $0.144\pm 0.002$ | $1.827\pm 0.004$
$p=0.40$ | $\lambda_{c}=0.241(0)$ | $0.883\pm 0.013$ | $0.208\pm 0.005$ | $1.856\pm 0.005$
$p=0.50$ | $\lambda_{c}=0.226(7)$ | $0.743\pm 0.009$ | $0.076\pm 0.002$ | $1.761\pm 0.004$
$p=0.60$ | $\lambda_{c}=0.213(3)$ | $0.775\pm 0.018$ | $0.100\pm 0.004$ | $1.780\pm 0.005$
$p=0.70$ | $\lambda_{c}=0.196(7)$ | $0.790\pm 0.016$ | $0.124\pm 0.004$ | $1.795\pm 0.005$
DP | $-$ | $3/4$ | $5/48$ | $43/24$

It is worth mentioning that MC time $t$ could be determined by incrementing
$t$ by $\delta t=1/n_{I}$, where $n_{I}$ represents the current number of
infected sites, each time an infected site is picked from the list. However,
in this context, we do not track time as our primary focus lies on static
quantities such as the fraction and mean cluster size of $R$ sites.

From the cluster size distribution, we have the fraction of recovery sites in
the finite (non-percolating) cluster with $s$ size

$P_{s}=s\frac{{n_{p}(s)}}{{n_{R}}},$ (3)

where $n_{R}$ is the total number of recovery sites and $n_{p}(s)$ is the
number of clusters with $s$ recovery sites. Furthermore, the fraction of
recovery sites in the percolating cluster $P_{\infty}$ can be obtained by

$P_{\infty}=1-\frac{1}{{n_{R}}}\sum\limits_{s}{s{\kern 1.0pt}n_{p}(s)},$ (4)

such that the above summation excludes the percolation cluster. Now we can
define the order parameter from Eq. (4) as

$P=<P_{\infty}>,$ (5)

where $<x>$ means an average taken over different dynamic realizations. The
epidemic phase of the model is reached when $P\neq 0$, that is, when the
percolating cluster density is non-zero; while the endemic phase happens when
$P=0$. Other important quantities are the mean cluster size

$S=\frac{1}{{n_{R}}}\sum\limits_{s}{s^{2}{\kern 1.0pt}n_{p}(s)},$ (6)

which plays the rule of the susceptibility in classical percolation theory
[31, 32] when taking the average over different runs, i.e.,

$\chi=<S>,$ (7)

the overall mean cluster size

$S^{\prime}=\frac{1}{{n_{R}}}\sideset{}{{}^{\prime}}{\sum}_{s}{s^{2}{\kern
1.0pt}n_{p}(s)},$ (8)

and the mean quadratic cluster size

$M^{\prime}=\frac{1}{{n_{R}}}\sideset{}{{}^{\prime}}{\sum}_{s}{s^{3}{\kern
1.0pt}n_{p}(s)},$ (9)

where the primed summations above also include the percolating cluster. It is
worth remarking that the last two quantities $S^{\prime}$ and $M^{\prime}$
only make sense for finite lattice as in the asymptotic limit
$(N\rightarrow\infty)$, the percolating cluster size diverges.

Figure 1: Static quantities for the stochastic SIR model on square lattices
with first- and second-neighbor interactions by assuming a dilution rate of
$p=0.20$. Panels (a), (b), and (c) display the order parameter $P$, Binder
cumulant $U$, and susceptibility $\chi$ as a function of the recovery rate
$\lambda$, respectively. From the Binder cumulant crossing, we can estimate
the epidemic threshold at $\lambda_{c}=0.258(1)$ for this case. Each data
point represents averages over $10^{5}$ independent realizations.

Figure 2: Static quantities for the stochastic SIR model on square lattices
with first- and second-neighbor interactions by assuming a dilution rate of
$p=0.50$. Panels (a), (b), and (c) display the order parameter $P$, Binder
cumulant $U$, and susceptibility $\chi$ as a function of the recovery rate
$\lambda$, respectively. From the Binder cumulant crossing, we can estimate
the epidemic threshold at $\lambda_{c}=0.226(7)$ for this case. Each data
point represents averages over $10^{5}$ independent realizations. Figure 3:
Plot for the critical threshold point for all considered $p$ cases. The red
curve represents the non-linear curve fitting according to
$\lambda_{c}=\lambda_{0}-a\exp(b\,p)$, being $\lambda_{0}=0.325$, $a=0.05$ and
$b=1.32$.

Figure 4: Panels (a), (b), and (c) display the log-log plot of the quantities
$\phi$, $P$, and $\chi$, respectively, for several bond-dilution cases ($p$)
calculated at $\lambda_{c}$ as a function of the linear size of the system
$L$. Straight lines represent the best linear fit to the corresponding data.

In the critical region, the cluster size distribution should obey a power-law
scaling [33, 34, 17] as

$n_{p}(s)=s^{-\tau}F[s^{\alpha}\;(\lambda-\lambda_{c})],$ (10)

where $\lambda_{c}$ is the epidemic threshold and $F$ is a scaling function.
Likewise, the remaining quantities also adhere to scaling relations as
prescribed by classical percolation theory.

$\displaystyle P$ $\displaystyle=L^{-\beta/\nu\;}\tilde{P}(L^{1/\nu}{\kern
1.0pt}|\lambda-\lambda_{c}|),$ (11) $\displaystyle\chi$
$\displaystyle=L^{\gamma/\nu\;}\tilde{\chi}(L^{1/\nu}{\kern
1.0pt}|\lambda-\lambda_{c}|),$ (12) $\displaystyle<S^{\prime}>$
$\displaystyle=L^{\gamma/\nu\;}\tilde{S}^{\prime}(L^{1/\nu}{\kern
1.0pt}|\lambda-\lambda_{c}|),$ (13) $\displaystyle<M^{\prime}>$
$\displaystyle=L^{(\beta+2\gamma)/\nu\;}\tilde{M}^{\prime}(L^{1/\nu}{\kern
1.0pt}|\lambda-\lambda_{c}|).$ (14)

The reciprocal correlation-length exponents $1/\nu$ can be obtained by
calculating the modulus of the logarithmic derivative of $P$ at the critical
threshold point $\lambda_{C}$

$\phi\equiv\bigg{|}\frac{d}{{d\lambda}}\ln(P)\bigg{|}_{\lambda=\lambda_{c}},$
(15)

where the derivative of the function $f\equiv\ln(P)$ was evaluated numerically
by using a finite central difference scheme in the form

$\frac{{df}}{{d\lambda}}\simeq\frac{1}{{2\,\Delta\lambda}}\left({f(\lambda+\Delta\lambda)-f(\lambda-\Delta\lambda)}\right),$
(16)

which has an truncation error of the order of
$\mathscr{O}(\Delta\lambda^{2})$. The error function $\delta f^{\prime}$ of
$df/d\lambda$ was obtained via error propagation from the uncertainties in the
values of $f$ ($\delta f$), being expressed by

$\delta f^{\prime}=\frac{1}{{2\,\Delta\lambda}}\sqrt{(\delta
f(\lambda+\Delta\lambda))^{2}+(\delta f(\lambda-\Delta\lambda))^{2}}.$ (17)

In our computations $\Delta\lambda=2.0\times 10^{-3}$. Close to $\lambda_{c}$,
the quantity $\phi$ obeys a power-law scaling as

$\phi=L^{1/\nu}\tilde{\phi}(L^{1/\nu}|\lambda-\lambda_{c}|),$ (18)

where $\tilde{\phi}$ is also a scaling function. Moreover, we can define a
universal quantity $U$ in which the scaling dependencies cancel out by
combining Eqs. (11), (13), and (14) in the following way [17]

$U=P\frac{{<M^{\prime}>}}{{<S^{\prime}>^{2}}},$ (19)

being analogous to the Binder cumulant for ferromagnetic spin model [35, 36],
and obeying also a scaling relation

$U=\tilde{U}\;(L^{1/\nu}{\kern 1.0pt}|\lambda-\lambda_{c}|).$ (20)

The crossing point of the $U$ curves for different lattice sizes allow us to
estimate the epidemic threshold $\lambda_{c}$. Once $\lambda_{c}$ is
determined, conducting finite-size scaling (FSS) analyses on the observables
$P$, $\chi$, and $\phi$ using Eqs. (11), (12), and (18) respectively, enables
the calculation of the corresponding critical exponent ratios $\beta/\nu$,
$\gamma/\nu$, and $1/\nu$. These ratios serve as indicators of the leading
critical exponents $\beta$, $\gamma$, and $\nu$, which collectively
characterize the universality class of the model.

The bond-diluted lattices utilized in this study were derived from regular
square lattices with free boundary conditions. Initially, we begin with a
standard square lattice where nodes are interconnected with their first and
second nearest neighbors through both outgoing and incoming links.
Subsequently, with a probability of $p$, we disconnect a selected site from
its respective nearest neighbors, ensuring at least one connection is retained
to prevent isolated sites. This process is iterated for each site, resulting
in a new lattice characterized by a node density of $q=1-p$ with connections
to both first and second nearest neighbors. Such networks emulate a more
realistic population, showcasing a diverse array of interconnected
relationships among individuals. The special case arises when $p=0$,
corresponding to pure lattices where all nodes maintain connectivity with
their first and second neighbors. This scenario was explored in detail in Ref.
[37], and it was found that this case belongs to the same universality class
as that of 2D direct percolation. It is important to note that only cases in
which $p\leq 0.70$ have been considered in the present study, as increasing
the bond dilution level beyond that could severally impact the cardinality of
the generated lattices.

Spanning clusters emerge when the cluster of recovered individuals extends
across the entire lattice, connecting opposite edges. This occurs when the
cloud spans from one side to the other. We generated over $10^{5}$ spanning
clusters for each $\lambda$ value considered to ensure accurate averaging of
the relevant quantities. We deal with several lattice sizes, ranging from
$N=400$ up to $N=25600$.

##  4  Results and Discussion

In this section, we show our numerical results of the stochastic SIR model
coupled to bond-diluted lattices. In order to determine both the critical
region and the order of the phase transition in this model on these lattices,
we calculated the order parameter $P$, Binder cumulant $U$, and susceptibility
$\chi$ for several $p$ cases in a wide range of the parameter $\lambda$. These
quantities were averaged over $10^{5}$ different dynamic realizations of the
SIR model. After identifying the critical region, we employed FSS analysis to
obtain reliable estimates of the critical $\lambda$ and the leading critical
exponents. For each $p$ case we consider different lattice configurations upon
taking the disorder averages.

In Figs. 1, 1 and 1 are shown the order parameter, Binder cumulant, and
susceptibility as a function of the recovery rate $\lambda$ for the case
$p=0.20$, respectively. While Figs. 2, 2 and 2 show the same quantities for
the case $p=0.50$. As one can see from Figs. 1 and 2, a typical second-order
phase transition take place for these cases in which a typical sigmoid-shaped
curve occurs. Similarly, the same conclusion can be drawn for the remaining
treated cases ($p\leq 0.70$). From Binder cumulant crossings, we can estimate
the corresponding epidemic thresholds for each case. The critical thresholds
were estimated with four significant figures. For the case $p=0.20$, we
obtained $\lambda_{c}=0.176(6)$, while for the case $p=0.50$, we got
$\lambda_{c}=0.226(7)$. Fig. 3 displays the critical threshold point for all
considered cases, including cases with high dilution degrees, namely $p=0.80$
and $p=0.90$. The red curve represents the best non-linear fit to the
corresponding data. One can observe from Fig. 3 that the critical point for
dilution disorder monotonically diminishes with increasing $p$ value.

By taking the slope of the log-log plot of the quantity $\phi$ versus the
linear size of the system $L$, we can estimate $1/\nu$ for each case. Fig. 4
shows the best linear fit to Eq. (18) for several cases. Error bars were
estimated by using Eq. (17). We obtained $1/\nu=0.667\pm 0.022$ for $p=0.10$,
$1/\nu=0.691\pm 0.019$ for $p=0.20$, $1/\nu=0.788\pm 0.011$ for $p=0.30$,
$1/\nu=0.883\pm 0.013$ for $p=0.40$, $1/\nu=0.743\pm 0.009$ for $p=0.50$,
$1/\nu=0.775\pm 0.018$ for $p=0.60$, and $1/\nu=0.790\pm 0.016$ for $p=0.70$.
As can be seen, the values of $1/\nu$ for the cases $p=0.10$, $p=0.20$,
$p=0.30$, and $p=0.40$ are several standard deviations away from the exact
critical exponent ratio $1/\nu=3/4$ of 2D dynamical percolation. Meanwhile,
for the remaining cases, the estimated values of $1/\nu$ are within a few
standard deviations of that exact value.

Similarly, a finite-size scaling analysis of the magnitudes of the order
parameter $P$ and the susceptibility $\chi$ at $\lambda_{c}$ by using Eqs.
(11) and (12) yield, respectively, the critical exponent ratios $\beta/\nu$
and $\gamma/\nu$ for all cases. Figs. 4 and 4 show the log-log plot of $P$ and
$\chi$ (both calculated at $\lambda_{c}$) against $L$, respectively. The
straight lines in those figures are the best linear fit to Eqs. (11) and (12),
respectively. The estimates for $\beta/\nu$ and $\gamma/\nu$ for all cases are
shown in Table 1.

Most of these values deviate significantly, by several standard deviations,
from the exact critical exponent ratios of 2D dynamical percolation,
specifically $\beta/\nu=5/48$ and $\gamma/\nu=43/24$. With the exception of
the case where $p=0.60$, all other estimates strongly suggest that the SIR
model on a bond-diluted lattice, in general, does not belong to the same
universality class as 2D dynamic percolation. Estimates of the critical
exponents ratios and critical threshold $\lambda_{c}$ for all cases are
summarized and compared with the corresponding exact values form 2D dynamics
percolation in Table 1.

##  5  Conclusions

We conducted Monte Carlo simulations of the stochastic SIR model on bond-
diluted square lattices to investigate the critical behavior exhibited by
these systems. We estimated both the critical threshold and the leading
critical exponent ratios for various cases. Our numerical analysis has
revealed that quenched bond dilution disorder significantly influences the
critical exponents of the model for most of the considered $p$ values. The
simulation results clearly suggest that the present model belongs to a new
universality class different from that of 2D dynamical percolation, depending
on the chosen $p$ value. Additionally, we observed that the critical threshold
in the stochastic SIR model monotonically decreases with increasing bond-
dilution disorder.

This study has broad relevance to many areas, not only to the dynamics of
infectious disease transmission but also to general diffusion processes, the
propagation of damage in random networks, and provides insights into how
perturbations in connectivity affect the overall behavior, resilience, and
robustness of complex systems. Additionally, the findings outlined in this
paper may contribute to a deeper understanding of how topological
irregularities can affect the critical properties of interconnected complex
systems.

##  6 Conflict of Interest Statement

The authors have no competing interests to declare that are relevant to the
content of this article.

##  7  References

## References

  * Kermack and Mckendrick [1927] W. . Kermack, A. G. Mckendrick,  A contribution to the mathematical theory of epidemics,  Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 115 (1927) 700–721. doi:10.1098/RSPA.1927.0118.
  * Kim et al. [2020] Y. Kim, A. V. Barber, S. Lee,  Modeling influenza transmission dynamics with media coverage data of the 2009 h1n1 outbreak in korea,  PLOS ONE 15 (2020) e0232580. doi:10.1371/JOURNAL.PONE.0232580.
  * Kassem and Ndam [2010] T. Kassem, J. Ndam,  A stochastic modeling of recurrent measles epidemic,  Science World Journal 3 (2010) 29–32. doi:10.4314/swj.v3i4.51824.
  * Tilahun et al. [2020] G. T. Tilahun, S. Demie, A. Eyob,  Stochastic model of measles transmission dynamics with double dose vaccination,  Infectious Disease Modelling 5 (2020) 478–494. doi:10.1016/J.IDM.2020.06.003.
  * Sun et al. [2023] X. Sun, C. Dai, K. Wang, Y. Liu, X. Jin, C. Wang, Y. Yin, Z. Ding, Z. Lu, W. Wang, Z. Wang, F. Tang, K. Wang, Z. Peng,  A dynamic compartmental model to explore the optimal strategy of varicella vaccination: An epidemiological study in jiangsu province, china,  Tropical Medicine and Infectious Disease 8 (2023) 17. doi:10.3390/TROPICALMED8010017/S1.
  * Huang et al. [2011] Z. Huang, Q. Yang, J. Cao,  Complex dynamics in a stochastic internal hiv model,  Chaos, Solitons & Fractals 44 (2011) 954–963. doi:10.1016/J.CHAOS.2011.07.017.
  * Wang et al. [2017] Y. Wang, D. Jiang, T. Hayat, B. Ahmad,  A stochastic hiv infection model with t-cell proliferation and ctl immune response,  Applied Mathematics and Computation 315 (2017) 477–493. doi:10.1016/j.amc.2017.07.062.
  * Atkeson et al. [2020] A. Atkeson, K. Kopecky, T. Zha,  Estimating and forecasting disease scenarios for covid-19 with an sir model,  NBER Working Paper No. 27335 (2020).
  * Yang et al. [2021] H. M. Yang, L. P. L. Junior, F. F. M. Castro, A. C. Yang,  Mathematical modeling of the transmission of sars-cov-2—evaluating the impact of isolation in são paulo state (brazil) and lockdown in spain associated with protective measures on the epidemic of covid-19,  PLOS ONE 16 (2021) e0252271. doi:10.1371/JOURNAL.PONE.0252271.
  * Gounane et al. [2021] S. Gounane, Y. Barkouch, A. Atlas, M. Bendahmane, F. Karami, D. Meskine,  An adaptive social distancing sir model for covid-19 disease spreading and forecasting,  Epidemiologic Methods 10 (2021) 20200044. doi:10.1515/em-2020-0044.
  * Wang et al. [2016] Z. Wang, C. T. Bauch, S. Bhattacharyya, A. d’Onofrio, P. Manfredi, M. Perc, N. Perra, M. Salathé, D. Zhao,  Statistical physics of vaccination,  Physics Reports 664 (2016) 1–113. doi:10.1016/j.physrep.2016.10.006.
  * Shao and Han [2022] Q. Shao, D. Han,  Epidemic spreading in metapopulation networks with heterogeneous mobility rates,  Applied Mathematics and Computation 412 (2022) 126559. doi:10.1016/J.AMC.2021.126559.
  * Helbing et al. [2015] D. Helbing, D. Brockmann, T. Chadefaux, K. Donnay, U. Blanke, O. Woolley-Meza, M. Moussaid, A. Johansson, J. Krause, S. Schutte, M. Perc,  Saving human lives: What complexity science and information systems can contribute,  Journal of Statistical Physics 158 (2015) 735–781. doi:10.1007/S10955-014-1024-9.
  * Grassberger [1983] P. Grassberger,  On the critical behavior of the general epidemic process and dynamical percolation,  Mathematical Biosciences 63 (1983) 157–172. doi:10.1016/0025-5564(82)90036-0.
  * Muñoz et al. [1999] M. A. Muñoz, R. Dickman, A. Vespignani, S. Zapperi,  Avalanche and spreading exponents in systems with absorbing states,  Physical Review E 59 (1999) 6175. doi:10.1103/PhysRevE.59.6175.
  * Tomé and Ziff [2010] T. Tomé, R. M. Ziff,  Critical behavior of the susceptible-infected-recovered model on a square lattice,  Physical Review E 82 (2010) 051921. doi:10.1103/PhysRevE.82.051921.
  * Souza et al. [2011] D. R. D. Souza, T. Tomé, R. M. Ziff,  A new scale-invariant ratio and finite-size scaling for the stochastic susceptible-infected-recovered model,  Journal of Statistical Mechanics: Theory and Experiment 2011 (2011) P03006. doi:10.1088/1742-5468/2011/03/P03006.
  * Pastor-Satorras et al. [2015] R. Pastor-Satorras, C. Castellano, P. V. Mieghem, A. Vespignani,  Epidemic processes in complex networks,  Reviews of Modern Physics 87 (2015) 925. doi:10.1103/RevModPhys.87.925.
  * Satulovsky and Tomé [1994] J. E. Satulovsky, T. Tomé,  Stochastic lattice gas model for a predator-prey system,  Physical Review E 49 (1994) 5073. doi:10.1103/PhysRevE.49.5073.
  * Antal et al. [2001] T. Antal, M. Droz, A. Lipowski, G. Ódor,  Critical behavior of a lattice prey-predator model,  Physical Review E 64 (2001) 036118. doi:10.1103/PhysRevE.64.036118.
  * Lugo and McKane [2008] C. A. Lugo, A. J. McKane,  Quasicycles in a spatial predator-prey model,  Physical Review E 78 (2008) 051911. doi:10.1103/PhysRevE.78.051911.
  * de Souza and Tomé [2010] D. R. de Souza, T. Tomé,  Stochastic lattice gas model describing the dynamics of the sirs epidemic process,  Physica A: Statistical Mechanics and its Applications 389 (2010) 1142–1150. doi:10.1016/j.physa.2009.10.039.
  * Newman and Ziff [2001] M. E. J. Newman, R. M. Ziff,  Fast monte carlo algorithm for site or bond percolation,  Physical Review E 64 (2001) 016706. doi:10.1103/PhysRevE.64.016706.
  * Ferraz and Herrmann [2008] C. H. A. Ferraz, H. J. Herrmann,  The strange man in random networks of automata,  Physica A: Statistical Mechanics and its Applications 387 (2008) 5689–5695. doi:10.1016/j.physa.2008.06.019.
  * Ferraz and Lima [2017] C. H. A. Ferraz, J. L. S. Lima,  Three-state potts model on non-local directed small-world lattices,  Physica A: Statistical Mechanics and its Applications 484 (2017) 488–498. doi:10.1016/j.physa.2017.05.016.
  * Ferraz and Lima [2018] C. H. A. Ferraz, J. L. S. Lima,  Kauffman cellular automata on quasicrystal topology,  Brazilian Journal of Physics 48 (2018) 531–538. doi:10.1007/s13538-018-0588-2.
  * Harris [1974] A. B. Harris,  Effect of random defects on the critical behaviour of ising models,  Journal of Physics C: Solid State Physics 7 (1974) 1671. doi:10.1088/0022-3719/7/9/009.
  * Barghathi and Vojta [2014] H. Barghathi, T. Vojta,  Phase transitions on random lattices: How random is topological disorder?,  Physical Review Letters 113 (2014) 120602. doi:10.1103/PHYSREVLETT.113.120602/FIGURES/4/MEDIUM.
  * Schrauth et al. [2018] M. Schrauth, J. S. E. Portela, F. Goth,  Violation of the harris-barghathi-vojta criterion,  Physical Review Letters 121 (2018) 100601. doi:10.1103/PhysRevLett.121.100601.
  * Sen [2011] P. Sen,  Probability distribution and sizes of spanning clusters at the percolation thresholds,  International Journal of Modern Physics C 8 (2011) 229–237. doi:10.1142/S0129183197000205.
  * Shante and Kirkpatrick [1971] V. K. Shante, S. Kirkpatrick,  An introduction to percolation theory,  Advances in Physics 20 (1971) 325–357. doi:10.1080/00018737100101261.
  * Stauffer and Aharony [2014] D. Stauffer, A. Aharony, Introduction to Percolation Theory, Taylor & Francis Group, 2014.
  * Stauffer [1979] D. Stauffer,  Scaling theory of percolation clusters,  Physics Reports 54 (1979) 1–74. doi:10.1016/0370-1573(79)90060-7.
  * Hoshen et al. [1979] J. Hoshen, D. Stauffer, G. H. Bishop, R. J. Harrison, G. D. Quinn,  Monte carlo experiments on cluster size distribution in percolation,  Journal of Physics A: Mathematical and General 12 (1979) 1285. doi:10.1088/0305-4470/12/8/022.
  * Selke and Shchur [2005] W. Selke, L. N. Shchur,  Critical binder cumulant in two-dimensional anisotropic ising models,  Journal of Physics A: Mathematical and General 38 (2005) L739. doi:10.1088/0305-4470/38/44/L03.
  * Ferraz [2015] C. H. A. Ferraz,  Critical behavior of the q = 3, 4 -potts model on quasiperiodic decagonal lattices,  Physica A: Statistical Mechanics and its Applications 440 (2015) 90–99. doi:10.1016/j.physa.2015.08.021.
  * Ferraz [2021] C. H. A. Ferraz,  A non-absorbing sir stochastic lattice gas model on hybrid lattices,  Physics Letters A 424 (2021) 127871. doi:10.1016/j.physleta.2021.127871.