# Efficient Modeling of Quasi-Periodic Data with Seasonal Gaussian Process

Ziang Zhang
Department of Statistical Sciences, University of Toronto
and
Patrick Brown
Department of Statistical Sciences, University of Toronto
Centre for Global Health Research, St Michael’s Hospital
and
Jamie Stafford
Department of Statistical Sciences, University of Toronto

###### Abstract

Quasi-periodicity refers to a pattern in a function where it appears periodic
but has evolving amplitudes over time. This is often the case in practical
settings such as the modeling of case counts of infectious disease or the
carbon dioxide (CO2) concentration over time. In this paper, we introduce a
class of Gaussian processes, called seasonal Gaussian Processes (sGP), for
model-based inference of such quasi-periodic behavior. We illustrate that the
exact sGP can be efficiently fit within $O(n)$ time using its state space
representation for equally spaced locations. However, for large datasets with
irregular spacing, the exact approach becomes computationally inefficient and
unstable. To address this, we develop a continuous finite dimensional
approximation for sGP using the seasonal B-spline (sB-spline) basis
constructed by damping B-splines with sinusoidal functions. We prove that the
proposed approximation converges in distribution to the true sGP as the number
of basis functions increases, and show its superior approximation quality
through numerical studies. We also provide a unified and interpretable way to
define priors for the sGP, based on the notion of predictive standard
deviation (PSD). Finally, we implement the proposed inference method on
several real data examples to illustrate its practical usage.

Keywords: Hierarchical Model, Bayesian Methods, Statistical Computing, Prior
Elicitation

##  Introduction

In astronomical studies the brightness of Sun-like stars tends to exhibit
seasonal variation with amplitudes evolving overtime owing to the rotational
modulation of magnetically active regions (Roberts et al.,, 2013). This kind
of quasi-periodic behavior can also be observed in the weekly mortality counts
due to influenza, and the carbon dioxide concentration (CO2) readings at the
Mauna Loa observatory (Rasmussen,, 2003). Therefore, inferring an unknown
function $g$ with quasi-periodic behavior is often a crucial task in a variety
of contexts.

In this paper, we proposed a class of Gaussian processes ($\mathcal{GP}$),
which we call seasonal Gaussian Processes (sGP), to make inference of a quasi-
periodic function $g$. The sGP has a natural construction through an ordinary
stochastic differential equation (SDE) which provides a clear interpretation
of the degree of quasi-periodicity in the inferred function. Additionally,
once written into a state space representation with its first derivative, the
sGP has a Markov property that simplifies the computation for inference at
equally spaced locations.

However, for large datasets with irregular spacing, the exact approach can
become computationally inefficient and numerically unstable. To address this,
we develop a novel finite dimensional continuous approximation to the sGP with
the seasonal B-spline basis (sB-spline). Motivated by the covariance structure
of the sGP, sB-splines are constructed by damping B-splines with sinusoidal
functions. The proposed sB-spline approximation is then constructed using the
finite element method (FEM)(Brenner et al.,, 2008) in a strategy similar to
Lindgren and Rue, (2008), Lindgren et al., (2011) and Zhang et al., 2023b  for
approximating different classes of Gaussian processes. Improvements in
computational efficiency are realized by adopting a least squares approach so
the resulting approximation is fully determined by a set of basis weights with
a highly sparse precision matrix. This approximation is shown to converge
rapidly to the true sGP as the number of basis functions increases.

The remainder of the paper is organized as follows. In Section 2, we introduce
the sGP model and its SDE characterization. We then derive some useful
properties of the sGP, in particular a state-space representation which makes
the sGP ideal for inference. In Section 3.1, we develop our proposed sB-spline
approximation to the sGP using the FEM. We also discuss the computational and
theoretical properties that justify its practical advantages. In Section 4, we
derive the predictive standard deviation (PSD) of the sGP (Zhang et al.,
2023b, ), and use it to define prior for the standard deviation parameter of
the sGP in a unified and interpretable way. In Section 5, we demonstrate the
effectiveness and versatility of our proposed methodology for inferring quasi-
periodic functions. We first apply the approach to three simple examples as a
proof of concept, and then showcase its advantage by re-analyzing the
atmospheric Carbon Dioxide (CO2) concentrations data collected in Hawaii,
which was originally analyzed in Rust et al., (1979). Our inference uses the
posterior approximation algorithm in Stringer et al., (2022); Bilodeau et al.,
(2022). For a more detailed description of how to implement this algorithm
with $\mathcal{GP}$, one can refer to Section 3 in Zhang et al., 2023b .

##  Inference with the Seasonal Gaussian Process

We consider the following hierarchical model:

$\displaystyle Y_{i}|\bm{\eta}$
$\displaystyle\overset{ind}{\sim}\pi(Y_{i}|\bm{\eta},\kappa),$ (1)
$\displaystyle\eta_{i}$
$\displaystyle=\bm{v}_{i}^{T}\bm{\beta}+\sum_{l=1}^{L}g_{l}(x_{li}),$
$\displaystyle g_{l}:$ $\displaystyle\Omega_{l}\rightarrow\mathbb{R},\
g_{l}\overset{ind}{\sim}\mathcal{GP}(\mathcal{C}_{l}),\forall l\in[L].$

Here $\pi(Y_{i}|\bm{\eta},\kappa)$ is a twice continuously-differentiable
density with linear predictors $\bm{\eta}=[\eta_{1},...,\eta_{n}]^{T}$,
covariates $\bm{v}_{i},\ x_{li}$ and hyperparameter $\kappa$. Each unknown
function $g_{l}$ is assigned an independent Gaussian Process ($\mathcal{GP}$)
model with zero mean and covariance function $C_{l}$. The model 1 is termed
the Extended Latent Gaussian Model in Stringer et al., (2022), and an example
is the partial likelihood of Cox proportional hazards model where each $Y_{i}$
depends on multiple elements of $\bm{\eta}$ (Zhang et al., 2023a, ). For the
purpose of exposition, we consider the inference of a single unknown function
$g$ with the sGP in the remainder of this section as well as Section 3 and
Section 4.

###  Seasonal Gaussian Process: A Differential Equation Characterization

The seasonal Gaussian process (sGP) is a zero-mean $\mathcal{GP}$ defined for
functions that exhibit quasi-periodic behavior. It is characterized by two
parameters, the frequency parameter $\alpha\in\mathbb{R}^{+}$ and the standard
deviation (SD) parameter $\sigma\in\mathbb{R}^{+}$. If the function $g$ is
assigned the sGP model, which we denote as $g\sim\text{sGP}(\alpha,\sigma)$,
then $g$ satisfies the following stochastic differential equation (SDE):

$\displaystyle
Lg(x)=\bigg{[}\frac{d^{2}}{dx^{2}}+\alpha^{2}\bigg{]}g(x)=\sigma\xi(x),$ (2)

where $\xi(x)$ is the Gaussian white noise process. Note that the null space
of the SDE in Eq. 2 is:

$\displaystyle\text{Null}\\{L\\}=\text{span}\\{\cos(\alpha x),\sin(\alpha
x)\\},$ (3)

which is the space of sinusoidal functions with period $2\pi/\alpha$. For
identifiability, we define the sGP $g$ in Eq. 2 to satisfy the following
initial conditions

$\displaystyle g(0)=g^{\prime}(0)=0,$ (4)

without the loss of generality. More general initial conditions can be
accommodated by adding two global trigonometric functions $\\{v_{1}\cos(\alpha
x),v_{2}\sin(\alpha x)\\}$ to $g$, where initial conditions become
$g(0)=v_{1}$ and $g^{\prime}(0)=v_{2}\alpha$.

The sGP model has a direct interpretation in its parameters: $\alpha$
specifies the frequency and the periodicity of the function $g$, and $\sigma$
specifies the degree that $g$ could deviate from the sinusoidal space in Eq.
3. When $\sigma$ is larger, the sample path of sGP can deviate more from the
sinusoidal pattern, which allows for more quasi-periodic behavior in the
inferred $g$. Therefore, the quasi-periodic behavior of $g$ such as in the
sunspot variation could be accommodated and quantified by the sGP model.

###  Properties of the Seasonal Gaussian Process

The specific covariance function of the sGP model in Eq. 2 is given in
Proposition 1.

######  Proposition 1 (Covariance Function of the Seasonal Gaussian Process).

Let $g\sim\text{sGP}(\alpha,\sigma)$. Then $g$ has a covariance function:

$\displaystyle C(x_{1},x_{2})$
$\displaystyle=\bigg{(}\frac{\sigma}{\alpha}\bigg{)}^{2}\bigg{[}\frac{x_{1}}{2}\cos(\alpha(x_{2}-x_{1}))-\frac{\cos(\alpha
x_{2})\sin(\alpha x_{1})}{2\alpha}\bigg{]}$ (5)
$\displaystyle=\bigg{(}\frac{\sigma}{\alpha}\bigg{)}^{2}\bigg{[}\frac{\cos(\alpha
x_{2})x_{1}}{2}\cos(\alpha x_{1})+\left(\frac{\sin(\alpha
x_{2})x_{1}}{2}-\frac{\cos(\alpha x_{2})}{2\alpha}\right)\sin(\alpha
x_{1})\bigg{]},$

for any $x_{1},x_{2}\in\mathbb{R}^{+}$ such that $x_{1}\leq x_{2}$.

The detailed derivation of Proposition 1 is provided in Appendix A. Since the
covariance function in Eq. 5 has the form of damped sinusoidal waves, it is
clear that the sGP should be viewed as a model for quasi-periodic functions.
As an illustration, some simulated sample paths from the sGP models with
different $\alpha$ are shown in Fig. 1.

(a) $\alpha=\pi$: Covariance (b) $\alpha=2\pi$: Covariance (c) $\alpha=\pi$:
Samples (d) $\alpha=2\pi$: Samples

Figure 1: (a,b) display the covariance functions of two sGPs at $3$, and (c,d)
display five sample paths from the two sGPs. The frequency parameter $\alpha$
equals to $\pi$ in (a,c) and $2\pi$ in (b,d), and the SD parameter $\sigma=1$
in both sGPs. Both the covariance functions and the sample paths exhibit
quasi-periodic behavior with amplitudes varying overtime.

Another property that makes sGP convenient to use in practice is that when
augmented with its first derivative, the resulting augmented model has a state
space representation that satisfies the Markov property. This property is
summarized in Theorem 1.

######  Theorem 1 (State Space Representation of the sGP).

Consider $g\sim\text{sGP}(\alpha,\sigma)$, and let
$\bm{s}=\\{s_{1},...,s_{n}\\}\subset\mathbb{R}^{+}$ denotes a set of $n$
sorted locations and spacing $d_{1}=s_{1}$ and $d_{i}=s_{i}-s_{i-1}$ for
$i\in\\{2,..,n\\}$. Then the augmented vector
$\bm{g}_{aug}(s_{i})=[g(s_{i}),g^{\prime}(s_{i})]^{T}$ can be written as a
Markov model:

$\addcontentsline{lla}{section}{\numberline{\string\crtrefnumber{equ:state-
space}}{e}qu:state-
space}\bm{{g}}_{aug}(s_{i+1})=\bm{R}_{i+1}\bm{{g}}_{aug}(s_{i})+\bm{\epsilon}_{i+1},$
(6)

where $\bm{\epsilon}_{i}\overset{ind}{\sim}N(0,\bm{\Sigma}_{i})$. The $2\times
2$ matrices $\bm{R}_{i}$ and $\bm{\Sigma}_{i}=\bm{Q}_{i}^{-1}$ are
respectively defined as:

$\displaystyle\bm{R}_{i}=\begin{bmatrix}\cos(\alpha
d_{i})&\frac{1}{\alpha}\sin(\alpha d_{i})\\\ -\alpha\sin(\alpha
d_{i})&\cos(\alpha d_{i})\\\ \end{bmatrix},\
\bm{\Sigma}_{i}=\sigma^{2}\begin{bmatrix}\frac{1}{\alpha^{2}}\bigg{(}\frac{d_{i}}{2}-\frac{\sin(2\alpha
d_{i})}{4\alpha}\bigg{)}&\frac{\sin^{2}(\alpha d_{i})}{2\alpha^{2}}\\\
\frac{\sin^{2}(\alpha d_{i})}{2\alpha^{2}}&\frac{2\alpha d_{i}+\sin(2\alpha
d_{i})}{4\alpha}\\\ \end{bmatrix}.$ (7)

A proof of theorem 1 can be found in Appendix B. An implication of Theorem 1
is that the discretely observed sGP model can be fitted sequentially using a
filtering method by augmenting the state space with first derivatives. This is
particularly useful when the response variable $y$ is Gaussian, as
demonstrated by Ansley and Kohn, (1990) and Ansley et al., (1993). Another
implication of Theorem 1 is that unlike
$[g(s_{1}),...,g(s_{n})]^{T}\in\mathbb{R}^{n}$ which has a dense precision
matrix, the augmented vector
$[\bm{g}_{aug}(s_{1})^{T},..,\bm{g}_{aug}(s_{n})^{T}]^{T}\in\mathbb{R}^{2n}$
will be a GMRF (Rue and Held,, 2005) with a larger but much sparser precision
matrix $\bm{Q}_{aug}$ being:

$\displaystyle\bm{Q}_{aug}=\left[\begin{smallmatrix}\bm{Q}_{1}+\bm{A}_{1}&\bm{H}_{1}&&&&\\\
\bm{H}_{1}^{T}&\bm{Q}_{1}+\bm{A}_{2}&\bm{H}_{2}&&&\\\
&\bm{H}_{2}^{T}&\bm{Q}_{2}+\bm{A}_{3}&\bm{H}_{3}&&\\\
&&\bm{H}_{3}^{T}&\bm{Q}_{3}+\bm{A}_{4}&\bm{H}_{4}&\\\ &&&&\ddots&&&\\\
&&&&&\bm{Q}_{n-2}+\bm{A}_{n-1}&\bm{H}_{n-1}\\\
&&&&\ldots&\bm{H}_{n-1}^{T}&\bm{Q}_{n-1}\\\ \end{smallmatrix}\right],$ (8)

where $\bm{Q}_{i}=\bm{\Sigma}_{i}^{-1}$,
$\bm{A}_{i}=\bm{R}_{i}^{T}\bm{Q}_{i}\bm{R}_{i}$ and
$\bm{H}_{i}=-\bm{R}_{i}^{T}\bm{Q}_{i}$.

When all the locations are equally spaced, the Markov model (1) is reduced to
a first-order vector auto-regression (VAR) model which requires only one
coefficient matrix $\bm{R}$ and one error covariance $\bm{\Sigma}$ to be
computed in order to obtain the entire precision matrix $\bm{Q}_{aug}$.
Because of the sparse structure of $\bm{Q}_{aug}$, the entire augmented vector
$[\bm{g}_{aug}(s_{1})^{T},..,\bm{g}_{aug}(s_{n})^{T}]^{T}$ can then be
efficiently incorporated into the modern approximate Bayesian inference
methods of Stringer et al., (2022); Rue et al., (2009). For large and
irregularly spaced datasets, however, the approach using the state-space
representation becomes computationally challenging, in particular when the
correlation is strong. This is because of the need to compute and invert each
matrix $\bm{\Sigma}_{i}$ in Eq. 8. It is therefore often useful to consider a
continuous finite-dimensional approximation to the entire process, which we
introduce in the following section.

##  Finite Dimensional Approximation to SGP

In this section, we derive the proposed seasonal B-spline (sB-spline)
approximation as a Least Square approximation through the Finite Element
Method (FEM). Our results demonstrate that this approximation has both
desirable computational and theoretical properties. In Section 3.1, we provide
the context and a general form of the FEM approximation. For the purpose of
computation efficiency, we then adopt a least square approach to construct the
FEM approximation and choose the cubic B-spline functions as the basis in
Section 3.2. In Section 3.3, we improve the FEM approximation by damping the
basis functions with sinusoidal functions, which finally yields our proposed
sB-spline approximation.

###  Finite Element Method

Given $\Omega=[a,b]$ for some $a,b\in\mathbb{R}^{+}$, the idea of FEM is to
construct an approximation $\tilde{g}_{k}$ that has the form of:

$\displaystyle\tilde{g}_{k}(x)=\sum_{i=1}^{k}w_{i}\varphi_{i}(x),$ (9)

where $\bm{w}=(w_{1},...,w_{k})^{T}\in\mathbb{R}^{k}$ is a set of random
weights and $\mathbb{B}_{k}=\\{\varphi_{i}(x),i\in[k]\\}$ is a set of $k$
predetermined basis functions. Given a set of test functions
$\mathbb{T}_{k}=\\{\phi_{i}(x);i\in[k]\\}$, the distribution of the unknown
weight vector $\bm{w}$ is defined such that

$\displaystyle\langle L\tilde{g}_{k}(x),\phi_{i}(x)\rangle$
$\displaystyle\overset{d}{=}\sigma\langle Lg(x),\phi_{i}(x)\rangle$ (10)
$\displaystyle=\sigma\langle\xi(x),\phi_{i}(x)\rangle,$

for any test function $\phi_{i}\in\mathbb{T}_{k}$, where $\langle
f_{1}(x),f_{2}(x)\rangle:=\int_{\Omega}f_{1}(x)f_{2}(x)dx$ denotes inner-
product between any $f_{1},f_{2}\in L^{2}$. The requirement Eq. 10 can be
vectorized as

$\displaystyle\textbf{B}\bm{w}\overset{d}{=}N(\bm{0},\sigma^{2}\textbf{T}),$
(11)

where B and T are $k\times k$ matrices with entries:

$\displaystyle\textbf{B}_{ij}=\langle\phi_{i}(x),L\varphi_{j}(x)\rangle,\quad\textbf{T}_{ij}=\langle\phi_{i}(x),\phi_{j}(x)\rangle.$
(12)

When both B and T are non-singular, the requirement in Eq. 10 implies
$\bm{w}\sim N(\bm{0},\sigma^{2}\bm{\Sigma}_{\bm{w}})$ where the precision
matrix $\bm{\Sigma}^{-1}_{\bm{w}}=\textbf{B}^{T}\textbf{T}^{-1}\textbf{B}$.

###  Least Square Approximation with Cubic B-Spline Basis

Given a set of locations $\bm{s}$, the FEM approximation at these locations
can be written as $\tilde{\bm{g}}_{k}(\bm{s})=\bm{\Phi}\bm{w}$ where
$\bm{\Phi}$ is $n\times k$ design matrix with $ij$ element
$\bm{\Phi}_{ij}=\varphi_{j}(s_{i})$. The computational bottleneck of the above
FEM approximation therefore depends on both $\bm{\Sigma}^{-1}_{\bm{w}}$ and
$\bm{\Phi}$. We aim to find a FEM approximation such that both
$\bm{\Sigma}^{-1}_{\bm{w}}$ and $\bm{\Phi}$ can be computed efficiently and
have sparse structures.

To simplify the computation of the precision matrix
$\bm{\Sigma}^{-1}_{\bm{w}}$, we consider a least square approach to construct
the FEM approximation, which corresponds to setting the test function
$\phi_{i}(x)=L\varphi_{i}(x)$ for each $i\in[k]$. This strategy reduces the
precision matrix to
$\bm{\Sigma}^{-1}_{\bm{w}}=\textbf{T}^{T}\textbf{T}^{-1}\textbf{T}=\textbf{T}$,
and hence avoids the need to compute another matrix B or to invert the matrix
T (Zhang et al., 2023b, ).

It remains to choose a set of basis functions $\mathbb{B}_{k}$ such that the
corresponding matrices T and $\bm{\Phi}$ are sparse. Given any basis
$\mathbb{B}_{k}$, it can be shown that the precision matrix T can be computed
as:

$\displaystyle\textbf{T}=a^{4}\textbf{G}+\textbf{C}+a^{2}\textbf{M},$ (13)

with $\textbf{G}_{ij}=\langle\varphi_{i},\varphi_{j}\rangle$,
$\textbf{C}_{ij}=\langle\frac{d^{2}\varphi_{i}}{dx^{2}},\frac{d^{2}\varphi_{j}}{dx^{2}}\rangle$
and
$\textbf{M}_{ij}=\langle\varphi_{i},\frac{d^{2}\varphi_{j}}{dx^{2}}\rangle+\langle\frac{d^{2}\varphi_{i}}{dx^{2}},\varphi_{j}\rangle$.
The detailed derivation is provided in Appendix C. Given Eq. 13, it is natural
to consider using the cubic B-spline basis as $\mathbb{B}_{k}$ in the
approximation, which ensures the sparsity of T with $O(8k)$ nonzero elements,
and the sparsity of $\bm{\Phi}$ with $O(4k)$ nonzero elements (de Boor,,
2000). The cubic B-spline basis yields a computationally efficient
approximation to the true sGP. Furthermore, the following Theorem 2 guarantees
the convergence of the approximation to the true sGP process as the number of
basis $k$ increases.

######  Theorem 2 (Convergence of B-spline Approximation).

Let $\Omega=[a,b]$ where $a,b\in\mathbb{R}^{+}$ and let
$g\sim\text{sGP}(\alpha,\sigma)$. Assume $\mathbb{B}_{k}$ is a set of $k$
cubic B-splines constructed with equally spaced knots over $\Omega$, and
$\tilde{g}_{k}$ denotes the corresponding FEM approximation defined as in Eq.
9, then:

$\lim_{k\to\infty}\mathcal{C}_{k}(x_{1},x_{2})=\mathcal{C}(x_{1},x_{2}),$

for any $x_{1},x_{2}\in\Omega$, where
$\mathcal{C}(x_{1},x_{2})=\operatorname{Cov}[g(x_{1}),g(x_{2})]$,
$\mathcal{C}_{k}(x_{1},x_{2})=\operatorname{Cov}[\tilde{g}_{k}(x_{1}),\tilde{g}_{k}(x_{2})]$.

The proof of this theorem is provided in Appendix D.

###  Improving Approximation Accuracy with sB-Spline

Although Theorem 2 guarantees the limiting property of the cubic B-spline
approximation, the error from the cubic B-spline approximation could be non-
negligible when only a limited number of basis is used, which is a typical
situation in real applications. To see that, note the covariance function (5)
at a fixed location $x_{1}$ can be written as:

$\displaystyle C_{x_{1}}(x)=C(x_{1},x)=l_{x_{1}}(x)\cos(\alpha
x)+\tilde{l}_{x_{1}}(x)\sin(\alpha x),$ (14)

where $l_{x_{1}}(x)$ and $\tilde{l}_{x_{1}}(x)$ are piecewise linear functions
with a single knot at $x_{1}$. In order to accurately approximate the
distribution of the original sGP, the basis $\mathbb{B}_{k}$ needs to be
flexible enough to approximate well the covariance function $C_{x_{1}}(x)$ for
any $x_{1}\in\Omega$. The cubic B-spline basis can well approximate the
covariance function when $\alpha$ is small relative to the width of $\Omega$.
However, when $\alpha$ is large relative to the width of $\Omega$, the
covariance function will be too oscillating to be approximated well by the
cubic B-spline basis. See Fig. 2 for example.

(a) $\alpha=2\pi/5$

(b) $\alpha=2\pi$

Figure 2: Assessing the approximation accuracy with cubic B-spline (red) and
sB-spline (green) basis. In order to make the y-axis comparable, we focus on
the correlation function $\rho_{5}(x)$ instead of the covariance function. The
y-axis shows the maximum correlation approximation error
$|\rho_{5}(x)-\tilde{\rho}_{5}(x)|$ for $x\in[1,9]$, where
$\rho_{5}(x)=\text{Corr}[g(5),g(x)]$ and
$\tilde{\rho}_{5}(x)=\text{Corr}[\tilde{g}_{k}(5),\tilde{g}_{k}(x)]$ are
respectively the correlation functions of the true sGP and of the
approximation. The interval $\Omega=[0,10]$, $\sigma=1$ and $\alpha=2\pi$ or
$2\pi/5$.

Given $C_{x_{1}}(x)$ is itself a piece-wise polynomial damped by sinusoidal
functions, we increase the accuracy of our approximation by augmenting the
B-spline basis with B-splines that are also damped by sinusoidal functions. We
call this new basis the seasonal B-spline (sB-spline) basis. The detailed
construction of the sB-spline is given as the following. We first construct
$r\in\mathbb{Z}^{+}$ cubic B-spline basis $\\{b_{i}(x)\\}_{i=1}^{r}$ over
$\Omega$. Then, we augment the cubic B-spline basis with its damped versions
and get

$\displaystyle\mathbb{B}_{k}=\\{b_{i}(x)\\}_{i=1}^{r}\cup\\{b_{i}(x)\cos(\alpha
x)\\}_{i=1}^{r}\cup\\{b_{i}(x)\sin(\alpha x)\\}_{i=1}^{r},$ (15)

where $k=3r$. Because the sB-spline basis is directly obtained by damping
B-spline it inherits the compact support of the B-spline basis (de Boor,,
2000). Therefore each of $\textbf{G},\textbf{C}$ and M is still a banded
sparse matrix, and so will be the precision matrix $\Sigma_{\bm{w}}^{-1}$ and
the design matrix $\bm{\Phi}$.

While preserving the original computational efficiency from cubic B-spline
basis, the sB-spline approximation has significantly better approximation
accuracy, particular when $k$ is small and $\alpha$ is large. This can be seen
from Fig. 2; When $\alpha=2\pi$, the sB-spline approximation reduces the
correlation error to less than $0.2$ with $k=20$, but the original cubic
B-spline requires $k$ more than $80$ to achieve the same accuracy. The sB-
spline approximation also inherits the convergence result in Theorem 2, as
summarized in the corollary below:

######  Corollary 1 (Convergence of sB-spline Approximation).

Let $\Omega=[a,b]$ where $a,b\in\mathbb{R}^{+}$ and let $g\sim\text{sGP}$ with
parameters $\alpha,\sigma>0$. Assume $\mathbb{B}_{k}$ is a set of $k$ sB-
splines constructed with equally spaced knots over $\Omega$ as defined in Eq.
15, and $\tilde{g}_{k}$ denotes the corresponding FEM approximation defined as
in Eq. 9, then:

$\lim_{k\to\infty}\mathcal{C}_{k}(x_{1},x_{2})=\mathcal{C}(x_{1},x_{2}),$

for any $x_{1},x_{2}\in\Omega$, where
$\mathcal{C}(x_{1},x_{2})=\operatorname{Cov}[g(x_{1}),g(x_{2})]$,
$\mathcal{C}_{k}(x_{1},x_{2})=\operatorname{Cov}[\tilde{g}_{k}(x_{1}),\tilde{g}_{k}(x_{2})]$.

###### Proof.

This corollary directly follows from Theorem 2 with the observation that Eq.
15 contains the original cubic B-spline basis. ∎

##  Prior for the sGP

To perform Bayesian inference with $g\sim\text{sGP}(\alpha,\sigma)$, it is
crucial to define a prior on the standard deviation (SD) parameter $\sigma$ in
an interpretable manner. However, the effect and interpretation of $\sigma$ on
the fitted function vary significantly with the choice of $\alpha$, even
though it can always be viewed as the deviation size from the corresponding
sinusoidal space $\text{Null}\\{L\\}$. Therefore, instead of assigning prior
on the original SD $\sigma$, we propose to assign the prior on the predictive
standard deviation (PSD) (Zhang et al., 2023b, ), which is a one-to-one
transformation of $\sigma$ with a consistent interpretation across different
values of $\alpha$.

Given a once-differentiable function $g$ and a prediction unit $h>0$, the PSD
is defined as:

$\addcontentsline{lla}{section}{\numberline{\string\crtrefnumber{equ:PSD-
general}}{e}qu:PSD-
general}\text{SD}\bigg{[}g(x+h)\bigg{|}g(x),g^{(1)}(x)\bigg{]},$ (16)

which quantifies the uncertainty in predicting the function $g$ at $h$ units
ahead, using the current value and the value of its derivative. In practice,
the unit $h$ can be chosen based on the practical requirements of the
application. For the sGP, the PSD has a simple formula, as given in the
following corollary:

######  Corollary 2 (Predictive Standard Deviation of the sGP).

Let $g\sim\text{sGP}$ with parameters $\alpha,\sigma>0$, then:

$\displaystyle\sigma(h)$
$\displaystyle=\text{SD}\bigg{[}g(x+h)\bigg{|}g(x),g^{(1)}(x)\bigg{]}=\frac{\sigma}{\alpha}\sqrt{\bigg{(}\frac{h}{2}-\frac{\sin(2\alpha
h)}{4\alpha}\bigg{)}},$ (17)

for any $x,h\in\mathbb{R}^{+}$.

###### Proof.

This corollary is a direct consequence of the state-space representation in
Theorem 1. ∎

It is worth noting that the PSD $\sigma(h)$ of the sGP does not depend on the
current location $x$, as shown in Eq. 17. This feature distinguishes our
approach from that of Sorbye and Rue, (2014), which uses the marginal standard
deviation $\text{SD}[g(x)]$ and hence requires a set of reference locations to
be specified. Additionally, the PSD in our approach is computed using the true
sGP and is therefore invariant to the approximation choices, such as the
number and placement of the knots. Finally, we choose an exponential prior for
$\sigma(h)$ that encourages $g$ to behave like a sinusoidal function, which is
motivated by Simpson et al., (2017). The prior for the original $\sigma$ is
then recovered by scaling the exponential prior for $\sigma(h)$ with
$\alpha/\sqrt{(h/2-\sin(2\alpha h)/4\alpha)}$.

##  Examples

###  Canadian Mortality

In this example, we utilize the sGP to analyze the cause-specific mortality
counts in Ontario from Statistics Canada. These data were collected weekly
from Jan 9th 2010 to January 1st 2022. To illustrate the methodology, we
consider the mortality counts due to heart attack, influenza and accidental
injuries. We fit models using data prior to March 1st, 2020 and use the
remaining data to quantify the excess mortality of each cause during the
COVID-19 pandemic.

For all three datasets (heart attack, influenza, and accidental injuries), we
use the following model:

$\displaystyle y_{i}|\lambda_{i}$
$\displaystyle\sim\text{Poisson}(\lambda_{i}),$ (18)
$\displaystyle\log(\lambda_{i})$
$\displaystyle=\eta_{i}=\beta_{0}+\beta_{1}x_{i}+g(x_{i})+\xi_{i},$
$\displaystyle g$ $\displaystyle\sim\text{sGP}(a=2\pi,\sigma),\ \xi_{i}\sim
N(0,\sigma_{\xi}^{2}),$

where $y_{i}$ denotes the mortality count and $x_{i}$ denotes the time of
measurement in years since Jan 9th 2010. The observation-level random effect
$\xi_{i}$ is used to capture the overdispersion. We assign independent normal
priors with mean 0 and variance 1000 to the fixed effects ($\beta_{0}$ and
$\beta_{1}$) and the boundary conditions of the sGP. For the variance
parameters $\sigma$ and $\sigma_{\xi}$, we use exponential priors such that:

$\displaystyle\mathbb{P}(\sigma(1)>0.01)=0.5,\ $
$\displaystyle\mathbb{P}(\sigma_{\xi}>1)=0.5,$ (19)

where $\sigma(1)$ is the $1$-year PSD of the sGP $g$ defined in Section 4.

To perform inference on the $1$-year sGP, we use the augmented space method
described in Section 2.1. The results of the analysis for the three mortality
rates are presented in Fig. 3. The mortality due to influenza shows a much
larger deviation from yearly periodicity compared to heart attack mortality.
In particular, the peak mortality rate in 2015 is much higher than peaks in
the preceding years. Consequently, the prediction uncertainty of $g$ is
greater for influenza compared to for heart attack. On the other hand, the
effect of seasonality in injury data is relatively small compared to that in
influenza and heart attack data. The inference results for total excess
mortality during the COVID-19 period are presented in Fig. 4. The posteriors
indicate strong evidence of negative excess mortality in influenza and likely
positive excess mortality in accidental injury. However, for heart attack,
there is no evidence of excess mortality in either direction.

Figure 3: Inferential Results for the mortality rates over time in Section
5.1. Figures (a)-(c) display the posterior distributions of
$\exp[\beta_{0}+\beta_{1}x_{i}+g(x_{i})]$ and the death counts (points). The
posterior means are shown in the blue solid lines and the 95 percent posterior
intervals are shown in the red dashed lines.

Figure 4: Inferential Results for total excess mortality in Section 5.1.
Figures (a)-(c) show the posterior samples of the total excess mortality,
where each sample is the sum of the difference between a posterior sample of
the predicted death and the actual death counts.

###  Lynx Counts

In this section, we analyze the annual counts of trapped lynxes in Canada for
1821–1934 (Campbell and Walker,, 1977), with the following model:

$\displaystyle y_{i}|\lambda_{i}$
$\displaystyle\sim\text{Poisson}(\lambda_{i}),\
\log(\lambda_{i})=\eta_{i}=\beta_{0}+g_{c}(x_{i})+g_{c/2}(x_{i})+\xi_{i},$
(20) $\displaystyle g_{j}$
$\displaystyle\sim\text{sGP}\bigg{(}a=\frac{2\pi}{j},\sigma_{j}\bigg{)}\text{
for }j\in\bigg{\\{}c,\frac{c}{2}\bigg{\\}},\ \xi_{i}\sim
N(0,\sigma_{\xi}^{2}),$

where $y_{i}$ denotes the counts of trapped lynxes in Canada and $x_{i}$
denotes the year since $1821$. We assume the variation in lynxes counts is
driven by two independent sGPs, with $c$ years ($\alpha=2\pi/c$) and $c/2$
years cycle ($\alpha=4\pi/c$). We assign independent exponential priors such
that $\mathbb{P}(\sigma_{j}(50)>1)=0.01$ for each $j\in\\{c,c/2\\}$, and
$\mathbb{P}(\sigma_{\xi}>1)=0.01$. All the boundary conditions in the sGPs and
the intercept $\beta_{0}$ are assigned with independent normal prior
$N(0,1000)$. The year of periodicity $c$ is assumed unknown between $6$ years
to $12$ years, and we place a discrete uniform prior for $c$ from $6$ to $12$
years with $0.1$ year spacing.

We carry out the inference of the sGPs with the state-space approach described
in Section 2.1. The inferential results are summarized in Fig. 5. Based on
Fig. 5(a), it is most likely that the data exhibits a roughly 10.1 years
cyclic variation, which is close to the duration of the solar cycle. The two
sGPs well capture the quasi-periodic behaviors in the variation, as
illustrated in Fig. 5(b).

Figure 5: Inferential Results for the Lynx analysis in Section 5.2. Figure (a)
shows the posterior distribution of the unknown periodicity year $c$. Figure
(b) shows the posterior distribution of $g_{c}+g_{c/2}$: posterior mean shown
in the blue line, 95 credible intervals shown in the red lines, and
observations shown in circles.

###  Sunspot Data

In this example, we analyze the yearly mean sunspot data from the World Data
Center SILSO, Royal Observatory of Belgium, Brussels (SILSO World Data
Center,, 2022), with the following model:

$\displaystyle y_{i}$ $\displaystyle=g_{s}(x_{i})+g_{tr}(x_{i})+\epsilon_{i},$
(21) $\displaystyle g_{s}$
$\displaystyle\sim\text{sGP}\bigg{(}a=\frac{2\pi}{c},\sigma_{s}\bigg{)},\
g_{tr}\sim\text{IWP}_{3}(\sigma_{tr}),$ $\displaystyle\epsilon_{i}$
$\displaystyle\sim N(0,\sigma_{\epsilon}^{2}),$

where $y_{i}$ denotes the log-transformed yearly mean sunspot, and $x_{i}$
denotes the year of measurement since 1700.

We model the seasonal variation in yearly mean sunspots with an sGP model with
periodicity around $c$ years, denoted as $g_{s}$. The long-term trend $g_{tr}$
is modelled with a third-order Integrated Wiener process (IWP-$3$) with
smoothing parameter $\sigma_{tr}$. For variance parameters, we use independent
exponential priors such that $\mathbb{P}(\sigma_{c}(50)>1)=0.5$,
$\mathbb{P}(\sigma_{tr}(50)>5)=0.5$ and $\mathbb{P}(\sigma_{\epsilon}>1)=0.1$.
All the boundary conditions in the sGP are assigned with independent normal
prior $N(0,1000)$. The year of periodicity $c$ is assumed unknown between $8$
years to $13$ years, and we place a discrete uniform prior for $c$ from $8$ to
$13$ years with $0.1$ year spacing.

To facilitate the computation, we then approximate $g_{tr}$ using 100 O-spline
basis functions introduced in (Zhang et al., 2023b, ), and $g_{s}$ using 300
sB-spline basis functions as described in Section 3.1. All the basis functions
are constructed with equally spaced knots. The inferential results are
summarized in Fig. 6. Similar to the results in Section 5.2, the posterior
distribution of $c$ is concentrated between $9.5$ years to $11.5$ years, close
to the duration of the solar cycle.

Figure 6: Inferential Results for the sunspot analysis in Section 5.3. Figure
(a) shows the posterior distribution of the unknown periodicity year $c$.
Figure (b) shows the posterior distribution of $g_{s}+g_{tr}$: posterior mean
shown in the blue line, 95 credible intervals shown in the red lines, and
observations shown in circles.

###  CO2 Variation

In this section, we will study the atmospheric Carbon Dioxide (CO2)
concentrations data collected from an observatory in Hawaii from 1960 to 2021.
The observations were collected monthly before May 1974 and weekly afterward.

Previous analyses have suggested the presence of several periodic components
in addition to the obvious yearly cycle in the CO2 variation. Rust et al.,
(1979) analyzed the data from this observatory and hypothesized that the
besides the obvious yearly cyclic variation, two longer periodic components
also likely exist, due to the southern oscillation index (ENSO) and solar
activity cycles. Therefore, Rust et al., (1979) adopted a model with an
exponential function for the trend, one sinusoidal function and its first
harmonic component for the annual cycle, and two additional sinusoidal
functions for the ENSO and solar cycles. In a later study of the climate
oscillation, Scafetta, (2010) also identified a 10.4-year cycle induced by the
solar cycle, and an additional 9.1-year cycle induced by long lunar tidal
behaviors

Subsequently, both Rasmussen, (2003) and Solin and Särkkä, (2014) analyzed the
CO2 data with Gaussian process regressions. They assumed the CO2 concentration
consists of an increasing long-term trend, a pronounced yearly seasonal
variation, and some irregular noise. Rasmussen, (2003) suggested that the
seasonal variation in the CO2 variation is smooth but not exactly periodic,
and both Rasmussen, (2003) and Solin and Särkkä, (2014) chose to model the
quasi-periodic seasonal variation as one GP with covariance being the product
of the canonical periodic covariance and the Matern covariance with $\nu=3/2$
or $\infty$.

Based on these existing studies, we consider the following model (M1):

$\displaystyle y_{i}$
$\displaystyle=g_{tr}(x_{i})+g_{1}(x_{i})+g_{\frac{1}{2}}(x_{i})+g_{\frac{44}{12}}(x_{i})+g_{9.1}(x_{i})+g_{10.4}(x_{i})+e_{i},$
(22) $\displaystyle g_{s}$
$\displaystyle\sim\text{sGP}\bigg{(}a=\frac{2\pi}{s},\sigma_{s}\bigg{)},\text{
for }s\in\bigg{\\{}1,\frac{1}{2},\frac{44}{12},9.1,10.4\bigg{\\}}$
$\displaystyle g_{tr}$ $\displaystyle\sim\text{IWP-3}(\sigma_{tr}),\ e_{i}\sim
N(0,\sigma_{e}^{2}),$

where each $y_{i}$ denotes the CO2 concentration and $x_{i}$ denotes the time
of measurement since March 30, 1960, in years. The long-term trend $g_{tr}$ is
modelled using a third-order IWP with smoothing parameter $\sigma_{tr}$. The
seasonal variation is modelled as the sum of five separate sGPs: the first two
sGPs capture the one-year cycle and its first-order harmonic, and the third
sGP capture the 44-month ENSO cycles suggested in Rust et al., (1979), and the
fourth and fifth sGPs respectively capture the 9.1-year lunar cycle and
10.4-year solar cycle suggested in Scafetta, (2010). All the components in the
model are considered to be independent.

In contrast to Rasmussen, (2003) and Solin and Särkkä, (2014) where all the
hyperparameters are estimated through maximizing marginal likelihood, we
assign independent exponential priors to all the variance parameters such
that:

$\displaystyle\mathbb{P}(\sigma_{s}(10)>1)=0.5,\text{ for
}s\in\bigg{\\{}1,\frac{1}{2},\frac{44}{12},9.1,10.4\bigg{\\}}$ (23)
$\displaystyle\mathbb{P}(\sigma_{tr}(10)>30)=0.5,\
\mathbb{P}(\sigma_{e}>1)=0.5.$

All the boundary conditions in the sGP or IWP are assigned with independent
priors $N(0,1000)$. As a comparison, we also consider an exact periodic model
(M2):

$\displaystyle y_{i}$
$\displaystyle=g_{tr}(x_{i})+g_{1}(x_{i})+g_{\frac{1}{2}}(x_{i})+g_{\frac{44}{12}}(x_{i})+g_{9.1}(x_{i})+g_{10.4}(x_{i})+e_{i},$
(24) $\displaystyle g_{s}(x)$ $\displaystyle=v_{s_{1}}\cos\bigg{(}\frac{2\pi
x}{s}\bigg{)}+v_{s_{2}}\sin\bigg{(}\frac{2\pi x}{s}\bigg{)},\text{ for
}s\in\bigg{\\{}1,\frac{1}{2},\frac{44}{12},9.1,10.4\bigg{\\}}$ $\displaystyle
g_{tr}$ $\displaystyle\sim\text{IWP-3}(\sigma_{tr}),\ e_{i}\sim
N(0,\sigma_{e}^{2}),$

where all the sGP components in M1 are replaced with just their corresponding
sinusoidal functions, and all the priors are the same as in M1.

In order to improve computational efficiency, we adopt FEM approximations for
all the $\mathcal{GP}$s in M1 and M2, with basis functions constructed using
equally spaced knots. Specifically, we utilize the sB-splines approximation
with 90 basis functions for the sGP models, and the O-Spline approximation
(Zhang et al., 2023b, ) with 50 basis functions for the IWP model. Once the
models M1 and M2 are fitted, we use them to forecast CO2 variation up to the
year 2030. The inference and prediction results are summarized in Fig. 7, Fig.
8 and Table 1.

As shown in Fig. 7, the sGP-based model (M1) produces a much more stable
prediction than the exact-periodic model (M2), demonstrated by the width of
the posterior interval in the forecast. This is because, the value of
$\sigma_{tr}(10)$ in M2 is estimated to be much larger (Post Median: 100.6)
compared to in M1 (Post Median: 0.593), due to the failure of M1 to
accommodate the quasi-periodic behavior of the seasonal variation. This
difference can be further confirmed in Fig. 8, where the derivative of the
long-term trend has smooth sample paths in M1 but very wiggly sample paths in
M2.

| 1st Quantile | Median | 3rd Quantile
---|---|---|---
$\sigma_{tr}(10)$ | 0.4371 (102.05) | 0.6446 (120.81) | 0.9462 (139.49)
$\sigma_{1}(10)$ | 1.8179 (–) | 2.4542 (–) | 3.1709 (–)
$\sigma_{\frac{1}{2}}(10)$ | 1.8179 (–) | 2.4542 (–) | 3.1709 (–)
$\sigma_{\frac{44}{12}}(10)$ | 0.4004 (–) | 0.4514 (–) | 0.5019 (–)
$\sigma_{9.1}(10)$ | 0.004248 (–) | 0.01045 (–) | 0.02445 (–)
$\sigma_{10.4}(10)$ | 0.01090 (–) | 0.02035 (–) | 0.03335 (–)
$\sigma_{\epsilon}$ | 0.5870 (0.6099) | 0.5928 (0.6159) | 0.5988 (0.6220)
Table 1: Posterior summary of variance parameters for the CO2 example in
Section 5.4. Results from M2 are shown in parenthesis.

(a) M1: Overall mean CO2 concentration (b) M1: Seasonal variation (c) M2:
Overall mean CO2 concentration (d) M2: Seasonal variation

Figure 7: Inferential Results for the CO2 example in Section 5.4. (a,c)
display the posterior of the overall mean CO2 concentration over years (i.e.
$g_{tr}+\sum_{s}g_{s}$). (b,d) show the posterior of the overall seasonal
component in the CO2 concentration (i.e. $\sum_{s}g_{s}$). In (a)-(d), the
blue solid lines denote the posterior mean, and the red dashed lines denote
the 95 percent posterior credible interval. The prediction result from M1 is
more stable than that from M2, as shown by the width of the interval.

(a) M1: Long-term trend $g_{tr}$ (b) M1: Inference of $g^{\prime}_{tr}$ (c)
M2: Long-term trend $g_{tr}$ (d) M2: Inference of $g^{\prime}_{tr}$

Figure 8: Detailed inferential results for the CO2 example in Section 5.4.
(a,c) display the posterior distribution of the long-term trend; (b,d) display
the posterior distribution of the first derivative of the long-term trend. In
(a)-(d), the blue solid lines denote the posterior mean, and the red dashed
lines denote the 95 percent posterior credible interval. The derivative of
$g_{tr}$ is estimated to be much smoother in M1, whereas in M2 the result is
oscillating due to the unaccounted quasi-periodicity.

##  Discussion

In this paper, we introduced a class of Gaussian processes called seasonal
Gaussian processes (sGP) for making model-based inference on quasi-periodic
functions. We then proposed an accurate and computationally efficient sB-
spline approximation to the sGP, which can handle densely irregularly spaced
locations. The proposed approximation is shown to converge to the sGP as the
number of basis functions increases, both theoretically and numerically. We
implemented our proposed approaches in a flexible and scalable manner using
modern approximate Bayesian inference algorithms, as described in Stringer,
(2020), and demonstrated their practical usage in diverse real data
applications.

While our focus was on making fully Bayesian inference for the unknown quasi-
periodic function $g$, our proposed approaches could also be implemented in a
mixed effect model with parameters estimated through maximum likelihood
(Rasmussen,, 2003) or generalized cross-validation (Golub et al.,, 1979). An
interesting avenue for future work is to extend our approach to incorporate
adaptive modeling of the sGP, where the standard deviation parameter $\sigma$
varies over time. This could be achieved using a strategy similar to the one
described in Yue et al., (2014).

We have made the software that implements our proposed methods available on a
Github repository, which can be found at https://github.com/AgueroZZ/sGPfit.

SUPPLEMENTARY MATERIAL

R-package:

The R-package sGPfit that implements the proposed methodology in the paper can
be found at github.com/AgueroZZ/sGPfit.

Codes and Data:

All the data used in this paper are publicly available for download. The
mortality data in Section 5.1 is available at statcan.gc.ca. The lynx data in
Section 5.2 is available in the R package datasets (version 4.2.1). The
sunspot data in Section 5.3 is available at sidc.be/silso. The CO2 data in
Section 5.4 can be accessed from scrippsco2.ucsd.edu.

The codes to replicate each of the examples can be found at
github.com/AgueroZZ/sGPcode.

Supplement:

The online supplement includes an R-markdown tutorial that provides step-by-
step guidance for users to implement our proposed approach using the Lynx
dataset in R.

## References

  * Ansley and Kohn, (1990) Ansley, C. F. and Kohn, R. (1990).  Filtering and smoothing in state space models with partially diffuse initial conditions.  Journal of Time Series Analysis, 11(4):275–293.
  * Ansley et al., (1993) Ansley, C. F., Kohn, R., and Wong, C.-M. (1993).  Nonparametric spline regression with prior information.  Biometrika, 80(1):75–88.
  * Bilodeau et al., (2022) Bilodeau, B., Stringer, A., and Tang, Y. (2022).  Stochastic convergence rates and applications of adaptive quadrature in bayesian inference.  Journal of the American Statistical Association, pages 1–11.
  * Brenner et al., (2008) Brenner, S. C., Scott, L. R., and Scott, L. R. (2008).  The mathematical theory of finite element methods, volume 3.  Springer.
  * Campbell and Walker, (1977) Campbell, M. and Walker, A. (1977).  A survey of statistical work on the mackenzie river series of annual canadian lynx trappings for the years 1821-1934 and a new analysis.  Journal of the Royal Statistical Society: Series A (General), 140(4):411–431.
  * de Boor, (2000) de Boor, C. (2000).  A practical guide to splines.  Springer.
  * Golub et al., (1979) Golub, G. H., Heath, M., and Wahba, G. (1979).  Generalized cross-validation as a method for choosing a good ridge parameter.  Technometrics, 21(2):215–223.
  * Harvey et al., (1990) Harvey, A. C. et al. (1990).  Forecasting, structural time series models and the kalman filter.  Cambridge Books.
  * Lindgren and Rue, (2008) Lindgren, F. and Rue, H. (2008).  On the second-order random walk model for irregular locations.  Scandanavian Journal of Statistics, 35(4):691 – 700.
  * Lindgren et al., (2011) Lindgren, F., Rue, H., and Lindstrom, J. (2011).  An explicit link between gaussian fields and gaussian markov random fields: the stochastic partial differential equation approach.  Journal of the Royal Statistical Society, Series B (Statistical Methodology), 73(4):423–498.
  * Rasmussen, (2003) Rasmussen, C. E. (2003).  Gaussian processes in machine learning.  In Summer school on machine learning, pages 63–71. Springer.
  * Roberts et al., (2013) Roberts, S., Osborne, M., Ebden, M., Reece, S., Gibson, N., and Aigrain, S. (2013).  Gaussian processes for time-series modelling.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 371(1984):20110550.
  * Rue and Held, (2005) Rue, H. and Held, L. (2005).  Gaussian Markov Random Fields: Theory and Applications.  Chapman and Hall/CRC Press.
  * Rue et al., (2009) Rue, H., Martino, S., and Chopin, N. (2009).  Approximate bayesian inference for latent gaussian models by using integrated nested laplace approximations.  Journal of the Royal Statistical Society. Series B (Statistical Methodology), 71(2):319 – 392.
  * Rust et al., (1979) Rust, B. W., Rotty, R. M., and Marland, G. (1979).  Inferences drawn from atmospheric co2 data.  Journal of Geophysical Research: Oceans, 84(C6):3115–3122.
  * Särkkä and Solin, (2019) Särkkä, S. and Solin, A. (2019).  Applied stochastic differential equations, volume 10.  Cambridge University Press.
  * Scafetta, (2010) Scafetta, N. (2010).  Empirical evidence for a celestial origin of the climate oscillations and its implications.  Journal of Atmospheric and Solar-Terrestrial Physics, 72(13):951–970.
  * Schultz, (1969) Schultz, M. H. (1969).  Approximation theory of multivariate spline functions in sobolev spaces.  SIAM Journal on Numerical Analysis, 6(4):570–582.
  * SILSO World Data Center, (2022) SILSO World Data Center (1700-2022).  The International Sunspot Number.  International Sunspot Number Monthly Bulletin and online catalogue.
  * Simpson et al., (2017) Simpson, D., Rue, H., Martins, T. G., Riebler, A., and Sorbye, S. H. (2017).  Penalising model component complexity: A principled, practical approach to constructing priors.  Statistical Science, 32(1):1 – 28.
  * Solin and Särkkä, (2014) Solin, A. and Särkkä, S. (2014).  Explicit link between periodic covariance functions and state space models.  In Artificial Intelligence and Statistics, pages 904–912. PMLR.
  * Sorbye and Rue, (2014) Sorbye, S. H. and Rue, H. (2014).  Scaling intrinsic gaussian markov random field priors in spatial modelling.  Spatial Statistics, 8:39–51.
  * Stringer, (2020) Stringer, A. (2020).  Implementing approximate bayesian inference using adaptive quadrature: the aghq package.  arXiv:2101.04468 [stat.CO].
  * Stringer et al., (2022) Stringer, A., Brown, P., and Stafford, J. (2022).  Fast, scalable approximations to posterior distributions in extended latent gaussian models.  Journal of Computational and Graphical Statistics, pages 1–36.
  * Yue et al., (2014) Yue, Y. R., Simpson, D., Lindgren, F., and Rue, H. (2014).  Bayesian adaptive smoothing splines using stochastic differential equations.
  * (26) Zhang, Z., Stringer, A., Brown, P., and Stafford, J. (2023a).  Bayesian inference for cox proportional hazard models with partial likelihoods, nonlinear covariate effects and correlated observations.  Statistical Methods in Medical Research, 32(1):165–180.
  * (27) Zhang, Z., Stringer, A., Brown, P., and Stafford, J. (2023b).  Model-based smoothing with integrated wiener processes and overlapping splines.  arXiv preprint arXiv:2302.02053.

## Appendix A: Derivation of the sGP covariance

######  Proposition 1 (Covariance Function of the Seasonal Gaussian Process).

Let $g\sim\text{sGP}(\alpha,\sigma)$. Then $g$ has a covariance function:

$\displaystyle C(x_{1},x_{2})$
$\displaystyle=\bigg{(}\frac{\sigma}{\alpha}\bigg{)}^{2}\bigg{[}\frac{x_{1}}{2}\cos(\alpha(x_{2}-x_{1}))-\frac{\cos(\alpha
x_{2})\sin(\alpha x_{1})}{2\alpha}\bigg{]}$ (25)
$\displaystyle=\bigg{(}\frac{\sigma}{\alpha}\bigg{)}^{2}\bigg{[}\frac{\cos(\alpha
x_{2})x_{1}}{2}\cos(\alpha x_{1})+\left(\frac{\sin(\alpha
x_{2})x_{1}}{2}-\frac{\cos(\alpha x_{2})}{2\alpha}\right)\sin(\alpha
x_{1})\bigg{]},$

for any $x_{1},x_{2}\in\mathbb{R}^{+}$ such that $x_{1}\leq x_{2}$.

###### Proof.

It is obvious that the differential operator $L$ is linear. Define
$\bm{g}_{aug}(x)=(g(x),g^{\prime}(x))^{T}$ and therefore
$\bm{g^{\prime}}_{aug}(x)=(g^{\prime}(x),g^{\prime\prime}(x))^{T}$, then the
SDE can be rewritten in the vector form:

$\displaystyle\bm{g^{\prime}}_{aug}=\bm{F}\bm{g}_{aug}+\bm{L}W,$ (26)

where $\bm{F}=\begin{bmatrix}0&1\\\ -\alpha^{2}&0\\\ \end{bmatrix}$ and
$\bm{L}=\begin{bmatrix}0\\\ \sigma\\\ \end{bmatrix}$.

Using the result from Särkkä and Solin, (2019) (section 4.3), the solution of
the linear SDE can be written as:

$\displaystyle\bm{g}_{aug}(x)$
$\displaystyle=\exp(\bm{F}x)\bm{g}_{aug}(0)+\int_{0}^{x}\exp(\bm{F}(x-\tau))\bm{L}W(\tau)d\tau$
(27) $\displaystyle=\int_{0}^{x}\exp(\bm{F}(x-\tau))\bm{L}W(\tau)d\tau,$

where $\exp(\bm{F}x)$ denotes the matrix exponential defined as
$\exp(\bm{F}x)=\sum_{k}\frac{\bm{F}^{k}x^{k}}{k!}$.

Note that $\bm{F}^{2k}=(-\alpha^{2})^{k}\bm{I}$ and
$\bm{F}^{2k}=(-\alpha^{2})^{k}\bm{F}$. With Taylor series, the first component
of $\bm{g}_{aug}(x)$ can be therefore written as:

$\displaystyle{g}(x)=\int_{0}^{x}\frac{\sigma}{\alpha}\sin(a(x-\tau))W(\tau)d\tau.$
(28)

Assume arbitrary $0<x_{1}\leq x_{2}$, the covariance function can be computed
for $g$ as:

$\displaystyle k(x_{1},x_{1})$
$\displaystyle=\int_{0}^{x_{1}}\frac{\sigma}{\alpha}\sin(\alpha(x_{1}-\tau))\frac{\sigma}{a}\sin(\alpha(x_{2}-\tau))d\tau$
(29)
$\displaystyle=\bigg{(}\frac{\sigma}{\alpha}\bigg{)}^{2}\bigg{[}\frac{x_{1}}{2}\cos(\alpha(x_{2}-x_{1}))-\frac{\cos(\alpha
x_{2})\sin(\alpha x_{1})}{2\alpha}\bigg{]},$

using properties of Gaussian white noise (Harvey et al.,, 1990).

∎

## Appendix B: Proof of the State-Space Representation

######  Theorem 1 (State Space Representation of the sGP).

Consider $g\sim\text{sGP}(\alpha,\sigma)$, and let
$\bm{s}=\\{s_{1},...,s_{n}\\}\subset\mathbb{R}^{+}$ denotes a set of $n$
sorted locations and spacing $d_{1}=s_{1}$ and $d_{i}=s_{i}-s_{i-1}$ for
$i\in\\{2,..,n\\}$. Then the augmented vector
$\bm{g}_{aug}(s_{i})=[g(s_{i}),g^{\prime}(s_{i})]^{T}$ can be written as a
Markov model:

$\bm{{g}}_{aug}(s_{i+1})=\bm{R}_{i+1}\bm{{g}}_{aug}(s_{i})+\bm{\epsilon}_{i+1},$
(30)

where $\bm{\epsilon}_{i}\overset{ind}{\sim}N(0,\bm{\Sigma}_{i})$. The $2\times
2$ matrices $\bm{R}_{i}$ and $\bm{\Sigma}_{i}=\bm{Q}_{i}^{-1}$ are
respectively defined as:

$\displaystyle\bm{R}_{i}=\begin{bmatrix}\cos(\alpha
d_{i})&\frac{1}{\alpha}\sin(\alpha d_{i})\\\ -\alpha\sin(\alpha
d_{i})&\cos(\alpha d_{i})\\\ \end{bmatrix},\
\bm{\Sigma}_{i}=\sigma^{2}\begin{bmatrix}\frac{1}{\alpha^{2}}\bigg{(}\frac{d_{i}}{2}-\frac{\sin(2\alpha
d_{i})}{4\alpha}\bigg{)}&\frac{\sin^{2}(\alpha d_{i})}{2\alpha^{2}}\\\
\frac{\sin^{2}(\alpha d_{i})}{2\alpha^{2}}&\frac{2\alpha d_{i}+\sin(2\alpha
d_{i})}{4\alpha}\\\ \end{bmatrix}.$ (31)

###### Proof.

To show the above Markov representation, note that the value of $g(s_{i+1})$
given $g(s_{i})$ can be written similarly as (Särkkä and Solin,, 2019):

$\bm{g}_{aug}(s_{i+1})=\exp(\bm{F}d_{i+1})\bm{g}_{aug}(s_{i})+\int_{s_{i}}^{s_{i+1}}\exp(\bm{F}(s_{i+1}-\tau))\bm{L}W(\tau)d\tau.$

Recall that $\bm{F}^{2k}=(-a^{2})^{k}\bm{I}$ and
$\bm{F}^{2k+1}=(-a^{2})^{k}\bm{F}$, then apply the Taylor series expansion for
both components in the integral above. It then can be rewritten as:

$\displaystyle\bm{g}_{aug}(s_{i+1})$
$\displaystyle=\exp(\bm{F}d_{i+1})\bm{g}_{aug}(s_{i})+\int_{s_{i}}^{s_{i+1}}\exp(\bm{F}(s_{i+1}-\tau))\bm{L}W(\tau)d\tau$
(32)
$\displaystyle=\bm{R}_{i+1}\bm{g}_{aug}(s_{i})+\int_{s_{i}}^{s_{i+1}}\begin{bmatrix}\frac{1}{a}\sin\big{(}a(s_{i+1}-\tau)\big{)}\\\
\cos\big{(}a(s_{i+1}-\tau)\big{)}\end{bmatrix}\sigma W(\tau)d\tau$
$\displaystyle:=\bm{R}_{i+1}\bm{g}_{aug}(s_{i})+\bm{\epsilon}_{i+1}.$

Note that since each $\bm{\epsilon}_{i+1}$ involves integration at disjoint
intervals, their independence follows from the property of Gaussian white
noise (Harvey et al.,, 1990). To check its covariance matrix
$\bm{\Sigma_{i+1}}$, note that:

$\displaystyle\bm{\Sigma_{i+1}}$
$\displaystyle=\sigma^{2}\leavevmode\resizebox{404.71352pt}{}{$\begin{bmatrix}\int_{s_{i}}^{s_{i+1}}\frac{1}{a^{2}}\sin^{2}\big{(}a(s_{i+1}-\tau)\big{)}d\tau&\frac{1}{a}\int_{s_{i}}^{s_{i+1}}\sin\big{(}a(s_{i+1}-\tau)\big{)}\cos\big{(}a(s_{i+1}-\tau)\big{)}d\tau\\\
\frac{1}{a}\int_{s_{i}}^{s_{i+1}}\sin\big{(}a(s_{i+1}-\tau)\big{)}\cos\big{(}a(s_{i+1}-\tau)\big{)}d\tau&\int_{s_{i}}^{s_{i+1}}\cos^{2}\big{(}a(s_{i+1}-\tau)\big{)}\end{bmatrix}$}$
(33)
$\displaystyle=\sigma^{2}\begin{bmatrix}\frac{1}{a^{2}}\bigg{(}\frac{d_{i+1}}{2}-\frac{\sin(2ad_{i+1})}{4a}\bigg{)}&\frac{\sin^{2}(ad_{i+1})}{2a^{2}}\\\
\frac{\sin^{2}(ad_{i+1})}{2a^{2}}&\frac{2ad_{i+1}+\sin(2ad_{i+1})}{4a}\\\
\end{bmatrix},$

which completes the proof.

∎

## Appendix C: Finite Element Method

The Finite Element Method (FEM) used to construct the finite-dimensional
approximation can be understood as the following procedures.

Given the (linear) stochastic differential equation (SDE) that defines the sGP
model:

$Lg(x)=\sigma\xi(x),$

where $L=a^{2}+\frac{d^{2}}{dx^{2}}$ is a linear differential operator and
$\xi(x)$ is the standard Gaussian white noise process. Let
$\Omega\subset\mathbb{R}^{+}$ denotes a bounded interval of interest. Let
$\mathbb{B}_{k}:=\\{\varphi_{i},i\in[k]\\}$ denote the set of $k$ pre-
specified basis functions, and let $\mathbb{T}_{q}:=\\{\phi_{i},i\in[q]\\}$
denote the set of $q$ pre-specified test functions. We consider finite
dimensional approximation with form
$\tilde{g}(.)=\sum_{i=1}^{k}w_{i}\varphi_{i}(.)$. The weights
$\bm{w}:=[w_{1},...,w_{k}]^{T}\in\mathbb{R}^{k}$ is a set of random weights to
be determined.

In our FEM construction, we used the sB-splines defined over $\Omega$ as the
basis functions, and chose the test functions by
$\mathbb{T}_{k}:=\\{\phi_{i}=L\varphi_{i},i\in[k]\\}$, which is called a least
squares approximation in Lindgren et al., (2011). The distribution of the
unknown weight vector can be found by fulfilling the weak formulation at the
test function spaces $\mathbb{T}_{k}$, such that

$\langle
L\tilde{g}(x),\phi_{i}(x)\rangle\overset{d}{=}\sigma\langle\xi(x),\phi_{i}(x)\rangle,$
(34)

for any test function $\phi_{i}\in\mathbb{T}_{k}$. This equation can also be
vectorized as:

$\langle L\tilde{g}(x),\phi_{i}(x)\rangle_{i=1}^{k}=H\bm{w},$

where the $ij$ component of the $k\times k$ $H$ matrix can be computed as
$H_{ij}=\langle L\varphi_{j}(x),L\varphi_{i}(x)\rangle_{i=1}^{k}.$

The inner product on the right $\langle\xi(x),\phi_{i}(x)\rangle_{i=1}^{k}$
will have Gaussian distribution with zero mean vector and covariance matrix
$H$ by properties of Gaussian white noise (Harvey et al.,, 1990). Therefore,
the basis coefficients $\bm{w}$ will be multivariate Gaussian with zero mean
and covariance $H^{-1}HH^{-1}=H^{-1}$. Each element of the matrix $H$ can be
written as:

$\displaystyle H_{ij}$ $\displaystyle=\langle
L\varphi_{j},L\varphi_{i}\rangle$ (35) $\displaystyle=\langle
a^{2}\varphi_{j}+\frac{d^{2}\varphi_{j}}{dx^{2}},a^{2}\varphi_{i}+\frac{d^{2}\varphi_{i}}{dx^{2}}\rangle$
$\displaystyle=a^{4}\langle\varphi_{j},\varphi_{i}\rangle+a^{2}\langle\frac{d^{2}\varphi_{j}}{dx^{2}},\varphi_{i}\rangle+a^{2}\langle\varphi_{j},\frac{d^{2}\varphi_{i}}{dx^{2}}\rangle+\langle\frac{d^{2}\varphi_{j}}{dx^{2}},\frac{d^{2}\varphi_{i}}{dx^{2}}\rangle,$

hence $H=a^{4}G+C+a^{2}M$ with $G_{ij}=\langle\varphi_{i},\varphi_{j}\rangle$,
$C_{ij}=\langle\frac{d^{2}\varphi_{i}}{dx^{2}},\frac{d^{2}\varphi_{j}}{dx^{2}}\rangle$
and
$M_{ij}=\langle\varphi_{i},\frac{d^{2}\varphi_{j}}{dx^{2}}\rangle+\langle\frac{d^{2}\varphi_{i}}{dx^{2}},\varphi_{j}\rangle$
for each element of the matrices.

## Appendix D: Proof of the Convergence Result

######  Theorem (Convergence of B-spline Approximation).

Let $\Omega=[a,b]$ where $a,b\in\mathbb{R}^{+}$ and let
$g\sim\text{sGP}(\alpha,\sigma)$. Assume $\mathbb{B}_{k}$ is a set of $k$
cubic B-splines constructed with equally spaced knots over $\Omega$, and
$\tilde{g}_{k}$ denotes the corresponding FEM approximation, then:

$\lim_{k\to\infty}\mathcal{C}_{k}(x_{1},x_{2})=\mathcal{C}(x_{1},x_{2}),$

for any $x_{1},x_{2}\in\Omega$, where
$\mathcal{C}(x_{1},x_{2})=\operatorname{Cov}[g(x_{1}),g(x_{2})]$,
$\mathcal{C}_{k}(x_{1},x_{2})=\operatorname{Cov}[\tilde{g}_{k}(x_{1}),\tilde{g}_{k}(x_{2})]$.

###### Proof.

The proof of this theorem starts with a similar strategy as in Lindgren et
al., (2011). Without the loss of generality, we assume the variance parameter
of the sGP $\sigma=1$ and $\Omega=[0,1]$. We denote the 3rd order Sobolev
space as
$H^{3}(\Omega)=\\{f\in\mathcal{L}^{2}(\Omega):D^{q}f\in\mathcal{L}^{2}(\Omega)\
\forall|q|\leq 3\\}$. We then define the constrained Sobolev space
$\mathcal{H}$ as:

$\mathcal{H}=\\{f\in H^{3}(\Omega):f(0)=f^{\prime}(0)=0\\}.$

Since $Lf=0$ implies $f\in\text{span}\\{\cos(\alpha x),\sin(\alpha x)\\}$, it
is clear that

$\displaystyle\left\langle f,h\right\rangle_{\mathcal{H}}:=\left\langle
Lf,Lh\right\rangle_{\Omega}=\int_{\Omega}Lf(x)Lh(x)dx$ (36)

defines an inner product for $f,h\in\mathcal{H}$.

Let $b>1$ and $\tilde{\Omega}=[0,b)\supset\Omega$, then for each $f\in
H^{3}(\Omega)$, there exists a zero-extension $f_{0}\in H^{3}(\tilde{\Omega})$
such that $f_{0}(x)=f(x)$ for $x\in\Omega$, and
$f_{0}(b)=f^{\prime}_{0}(b)=0$. Hence, we assume without the loss of
generality that for each $f\in\mathcal{H}$, $f(1)=f^{\prime}(1)=0$ from now.
This implies that the differential operator $D^{q}$ in $\mathcal{H}$ has
adjoint operator $(D^{q})^{*}=(-1)^{q}D^{q}$ for each $q\in\mathbb{Z}$.

Define $\mathcal{H}_{k}=\text{span}\\{\mathbb{B}_{k}\\}$. Note
$\mathcal{H}_{k}\subset\mathcal{H}$ by our construction of the B-spline basis.
This implies if $f(x)\in\mathcal{H}$, then there exists a projection
$\tilde{f}(x)=\sum_{i=1}^{k}w_{i}\varphi_{i}(x)\in\mathcal{H}_{k}$ which
satisfies:

$\displaystyle\left\langle
f-\tilde{f},\tilde{h}\right\rangle_{\mathcal{H}}=\left\langle
f,\tilde{h}\right\rangle_{\mathcal{H}}-\left\langle\tilde{f},\tilde{h}\right\rangle_{\mathcal{H}}=0,\
\forall\tilde{h}\in\mathcal{H}_{k}.$ (37)

To prove
$\lim_{k\to\infty}\mathcal{C}_{k}(x_{1},x_{2})=\mathcal{C}(x_{1},x_{2})$, note
that for any $f,h\in\mathcal{H}$,

$\displaystyle\operatorname{Cov}[\left\langle
g,f\right\rangle_{\mathcal{H}},\left\langle g,h\right\rangle_{\mathcal{H}}]$
$\displaystyle:=\operatorname{Cov}[\left\langle
Lg,Lf\right\rangle_{\Omega},\left\langle Lg,Lh\right\rangle_{\Omega}]$ (38)
$\displaystyle=\operatorname{Cov}[\left\langle\xi,Lf\right\rangle_{\Omega},\left\langle\xi,Lh\right\rangle_{\Omega}]$
$\displaystyle=\left\langle Lf,Lh\right\rangle_{\Omega},$
$\displaystyle=\left\langle f,h\right\rangle_{\mathcal{H}}.$

The second equality follows from the definition of the sGP, and the third
equality follows from the property of Gaussian white noise (Harvey et al.,,
1990).

Since the B-spline approximation $\tilde{g}_{k}\in\mathcal{H}_{k}$, we know

$\left\langle\tilde{g}_{k},f\right\rangle_{\mathcal{H}}=\left\langle\tilde{g}_{k},f-\tilde{f}+\tilde{f}\right\rangle_{\mathcal{H}}=\left\langle\tilde{g}_{k},\tilde{f}\right\rangle_{\mathcal{H}}+\left\langle\tilde{g}_{k},f-\tilde{f}\right\rangle_{\mathcal{H}}=\left\langle\tilde{g}_{k},\tilde{f}\right\rangle_{\mathcal{H}}.$

Using this result and the fact that the B-spline approximation $\tilde{g}_{k}$
is a least square solution, we have

$\displaystyle\operatorname{Cov}\bigg{[}\left\langle\tilde{g}_{k},f\right\rangle_{\mathcal{H}},\left\langle\tilde{g}_{k},h\right\rangle_{\mathcal{H}}\bigg{]}$
$\displaystyle=\operatorname{Cov}\bigg{[}\left\langle\tilde{g}_{k},\tilde{f}\right\rangle_{\mathcal{H}},\left\langle\tilde{g}_{k},\tilde{h}\right\rangle_{\mathcal{H}}\bigg{]}$
(39)
$\displaystyle=\operatorname{Cov}\bigg{[}\left\langle\xi,L\tilde{f}\right\rangle_{\Omega},\left\langle\xi,L\tilde{h}\right\rangle_{\Omega}\bigg{]}$
$\displaystyle=\left\langle\tilde{f},\tilde{h}\right\rangle_{\mathcal{H}}.$

Let $\mathcal{C}_{s}(x)=\mathcal{C}(s,x)$ denote the covariance function of
the sGP defined at any $s\in\Omega$. Based on the previous result in
Proposition 1, we know $\mathcal{C}_{s}(x)\in\mathcal{H}$ and
$L\mathcal{C}_{s}(x)=\frac{1}{\alpha}\sin[\alpha(s-x)^{+}]$ is the Green
function of $L$. Its projection into $\mathcal{H}_{k}$ is denoted as
$\tilde{\mathcal{C}}_{s}(x)$

######  Lemma 1.

Given the same setting in the main theorem

$\displaystyle\mathcal{C}(x_{1},x_{2})$
$\displaystyle=\operatorname{Cov}\big{[}\left\langle
g,{\mathcal{C}}_{x_{1}}\right\rangle_{\mathcal{H}},\left\langle
g,{\mathcal{C}}_{x_{2}}\right\rangle_{\mathcal{H}}\big{]}$ (40)
$\displaystyle\mathcal{C}_{k}(x_{1},x_{2})$
$\displaystyle=\operatorname{Cov}\bigg{[}\left\langle\tilde{g}_{k},\tilde{\mathcal{C}}_{x_{1}}\right\rangle_{\mathcal{H}},\left\langle\tilde{g}_{k},\tilde{\mathcal{C}}_{x_{2}}\right\rangle_{\mathcal{H}}\bigg{]}.$

###### Proof.

The first part directly follows from the proof in Proposition 1. The second
part can be proved using the fact that $L$ is self-adjoint and
$L\mathcal{C}_{x_{1}}(x)=\frac{1}{a}\sin[a(x_{1}-x)^{+}]$ is the Green
function, which implies
$\left\langle\varphi_{i},\tilde{\mathcal{C}}_{x_{1}}\right\rangle_{\mathcal{H}}=\varphi_{i}(x_{1})$
for each $i\in[k]$. The detailed proof proceeds as the following.

By construction of the B-spline approximation,

$\displaystyle\mathcal{C}_{k}(x_{1},x_{2})$
$\displaystyle=\text{Cov}\bigg{[}\sum_{i}w_{i}\varphi_{i}(x_{1}),\sum_{i}w_{i}\varphi_{i}(x_{2})\
\bigg{]}$ (41)
$\displaystyle=\text{Cov}\bigg{[}\bm{\Phi}(x_{1})^{T}\bm{w},\bm{\Phi}(x_{2})^{T}\bm{w}\bigg{]}$
$\displaystyle=\bm{\Phi}(x_{1})^{T}\Sigma_{\bm{w}}\bm{\Phi}(x_{2}),$

where $\bm{\Phi}(x)=[\varphi_{1}(x),...,\varphi_{k}(x)]^{T}$.

Since $\tilde{\mathcal{C}}_{x_{1}}(x)$ is the projection of
${\mathcal{C}}_{x_{1}}(x)$ to $\mathcal{H}_{k}$,
$\tilde{\mathcal{C}}_{x_{1}}(x)=\sum_{i}w_{x_{1},i}\varphi_{i}(x)$ for some
weights $\bm{w}_{x_{1}}=[w_{x_{1},1},...,w_{x_{1},k}]^{T}$. The same argument
can be used for $\tilde{\mathcal{C}}_{x_{2}}$. Therefore

$\displaystyle\text{Cov}\bigg{[}\left\langle\tilde{g}_{k},\tilde{\mathcal{C}}_{x_{1}}\right\rangle_{\mathcal{H}},\left\langle\tilde{g}_{k},\tilde{\mathcal{C}}_{x_{2}}\right\rangle_{\mathcal{H}}\bigg{]}$
$\displaystyle=\left\langle\tilde{\mathcal{C}}_{x_{1}},\tilde{\mathcal{C}}_{x_{2}}\right\rangle_{\mathcal{H}}$
(42)
$\displaystyle=\left\langle\bm{\Phi}(x)^{T}\bm{w}_{x_{1}},\bm{\Phi}(x)^{T}\bm{w}_{x_{2}}\right\rangle_{\mathcal{H}}$
$\displaystyle=\left\langle
L\bm{\Phi}(x)^{T}\bm{w}_{x_{1}},L\bm{\Phi}(x)^{T}\bm{w}_{x_{2}}\right\rangle_{\Omega}$
$\displaystyle=\bm{w}_{x_{1}}^{T}\Sigma^{-1}_{\bm{w}}\bm{w}_{x_{2}}.$

Solving
$\left\langle\tilde{\mathcal{C}}_{x_{1}},\varphi_{i}\right\rangle_{\mathcal{H}}=\left\langle\mathcal{C}_{x_{1}},\varphi_{i}\right\rangle_{\mathcal{H}}$
for each $i\in[k]$ yields that
$\bm{w}_{x_{1}}=\Sigma_{\bm{w}}\bm{\omega}_{x_{1}}$, where
$\bm{\omega}_{x_{1}}\in\mathbb{R}^{k}$ with $i$th element
${\omega}_{x_{1},i}=\left\langle\mathcal{C}_{x_{1}},\varphi_{i}\right\rangle_{\mathcal{H}}$.
Hence it only remains to show $\bm{\omega}_{x_{1}}=\bm{\Phi}(x_{1})$, which
holds because

$\displaystyle\omega_{x_{1},i}$
$\displaystyle=\left\langle\mathcal{C}_{x_{1}},\varphi_{i}\right\rangle_{\mathcal{H}}$
(43) $\displaystyle=\left\langle
L\mathcal{C}_{x_{1}},L\varphi_{i}\right\rangle_{\Omega}$
$\displaystyle=\left\langle
L^{*}L\mathcal{C}_{x_{1}},\varphi_{i}\right\rangle_{\Omega}$
$\displaystyle=\varphi_{i}(x_{1}),\ \forall i\in[k].$

The last equality holds because $L$ is a self-adjoint operator and
$L\mathcal{C}_{x_{1}}$ is the Green function. This lemma is hence proved. ∎

Using the above result and Lemma 1, it suffices to prove that

$\displaystyle\left\langle\tilde{\mathcal{C}}_{x_{1}},\tilde{\mathcal{C}}_{x_{2}}\right\rangle_{\mathcal{H}}\rightarrow\left\langle\mathcal{C}_{x_{1}},\mathcal{C}_{x_{2}}\right\rangle_{\mathcal{H}},$
(44)

as $k\rightarrow\infty$. For this step, we will use the following lemma on the
spline approximation:

######  Lemma 2.

Given the same setting in the main theorem, define the norm
$||f||_{\mathcal{H}}=\left\langle f,f\right\rangle_{\mathcal{H}}^{1/2}$ for
$f\in\mathcal{H}$ then

$\displaystyle||\mathcal{C}_{s}-\tilde{\mathcal{C}}_{s}||_{\mathcal{H}}$
$\displaystyle=O(1/k),$ (45)

for each $s\in\Omega$.

###### Proof.

The proof of this lemma mostly follows from the result in Schultz, (1969).
Given $\mathcal{H}_{k}$ is a spline space with degree $3$ and mesh size $1/k$
and $\mathcal{C}_{s}\in H^{3}(\Omega)$, by theorem 3.3 in Schultz, (1969) we
have

$\displaystyle||D^{q}(\mathcal{C}_{s}-\tilde{\mathcal{C}}_{s})||_{H^{2}(\Omega)}\leq
c_{q}\bigg{(}\frac{1}{k}\bigg{)}^{3-q}$ (46)

for each $0\leq q\leq 2$, where $c_{q}$ is a constant that depends on
$||\mathcal{C}_{s}||_{H^{3}(\Omega)}$ but does not depend on $k$. Note that

$\displaystyle||\mathcal{C}_{s}-\tilde{\mathcal{C}}_{s}||^{2}_{\mathcal{H}}$
$\displaystyle=\alpha^{4}||\mathcal{C}_{s}-\tilde{\mathcal{C}}_{s}||^{2}_{\mathcal{L}^{2}}+||D^{2}(\mathcal{C}_{s}-\tilde{\mathcal{C}}_{s})||^{2}_{\mathcal{L}^{2}}-2\alpha\left\langle(\mathcal{C}_{s}-\tilde{\mathcal{C}}_{s}),D^{2}(\mathcal{C}_{s}-\tilde{\mathcal{C}}_{s})\right\rangle_{\Omega}$
(47)
$\displaystyle=\alpha^{4}||\mathcal{C}_{s}-\tilde{\mathcal{C}}_{s}||^{2}_{\mathcal{L}^{2}}+||D^{2}(\mathcal{C}_{s}-\tilde{\mathcal{C}}_{s})||^{2}_{\mathcal{L}^{2}}+2\alpha||D(\mathcal{C}_{s}-\tilde{\mathcal{C}}_{s})||^{2}_{\mathcal{L}^{2}},$

where the second equality holds since $D^{*}=-D$. The lemma is proved since
$||.||_{\mathcal{H}}$ and $||.||_{H^{2}(\Omega)}$ are equivalent norm.

∎

Using Lemma 2, Eq. 44 can be proved with Cauchy Schwarz inequality and the
fact that the sequence $||\tilde{\mathcal{C}}_{x_{2}}||_{\mathcal{H}}$ is
bounded,

$\displaystyle\bigg{|}\left\langle\tilde{\mathcal{C}}_{x_{1}},\tilde{\mathcal{C}}_{x_{2}}\right\rangle_{\mathcal{H}}-\left\langle\mathcal{C}_{x_{1}},\mathcal{C}_{x_{2}}\right\rangle_{\mathcal{H}}\bigg{|}$
$\displaystyle=\bigg{|}\left\langle\tilde{\mathcal{C}}_{x_{1}}-\mathcal{C}_{x_{1}},\tilde{\mathcal{C}}_{x_{2}}\right\rangle_{\mathcal{H}}-\left\langle\mathcal{C}_{x_{1}},\mathcal{C}_{x_{2}}-\tilde{\mathcal{C}}_{x_{2}}\right\rangle_{\mathcal{H}}\bigg{|}$
(48)
$\displaystyle\leq\big{|}\big{|}\tilde{\mathcal{C}}_{x_{1}}-\mathcal{C}_{x_{1}}\big{|}\big{|}_{\mathcal{H}}\big{|}\big{|}\tilde{\mathcal{C}}_{x_{2}}\big{|}\big{|}_{\mathcal{H}}+\big{|}\big{|}\tilde{\mathcal{C}}_{x_{2}}-\mathcal{C}_{x_{2}}\big{|}\big{|}_{\mathcal{H}}\big{|}\big{|}\mathcal{C}_{x_{1}}\big{|}\big{|}_{\mathcal{H}}$
$\displaystyle=c/k,$

where $c$ is some constant independent of $k$. The theorem is hence proved. ∎


The result then follows from noting that $\mu\leq 1$ and that the welfare of
$M$ is the sum of the expected utilities and revenue. ∎

###  5.2 Simultaneous Composition

We now apply our generalized framework and study the impact of local
mechanisms’ properties on the equilibrium performance of a larger system.
Consider $n$ single-parameter agents seeking abstract service. Agents may
participate in one or more of $m$ separate mechanisms. Each of the mechanisms
are run simultaneously; each agent takes a profile of actions, one per
mechanism. Taking allocation levels in $[0,1]$ to be fractional levels of
service, we define an agent’s service level in the combined mechanism to be
their maximum service level across all mechanisms. Agents make the assigned
payments to all mechanisms. In this section we show that under these
assumptions, if each individual mechanisms has weak competitive efficiency
$\mu$, then so too does the aggregate mechanism induced for agents by
simultaneous participation as described above. We refer to the aggregate
mechanism as the simultaneous composition of the individual component
mechanisms.

Formally, a simultaneous composition of mechanisms consists of $m$ separate
feasibility environments $\mathcal{X}^{1},\ldots,\mathcal{X}^{m}$, one per
component mechanism. Each component mechanism $M^{j}$ is comprised of a bid
allocation rule $\mathbf{\tilde{x}}^{j}$ and a bid payment rule
$\mathbf{\tilde{p}}^{j}$, mapping a profile of actions $\mathbf{a}^{j}$ to an
allocation in $\mathcal{X}^{j}$ and a nonnegative payment, respectively. We
assume each mechanism has a withdraw action $\bot$ which guarantees zero
allocation and payments in $M^{j}$. Define the simultaneous composition of
mechanisms $M^{1},\ldots,M^{m}$ in the following way:

######  Definition 70.

Let mechanisms $M^{1},\ldots,M^{m}$ have bid allocation and bid payment rules
$(\mathbf{\tilde{x}}^{j},\mathbf{\tilde{p}}^{j})$ and individual action spaces
spaces $A_{i}^{1},\ldots,A_{i}^{m}$ for each agent $i$. The _simultaneous
composition_ of $M^{1},\ldots,M^{m}$ is defined to have:

  * •

Action space $\prod_{j}A_{i}^{j}$ for each agent. That is, each agent
participates in the global mechanism by participating in each composed
mechanism individually. Given a profile of actions $\mathbf{a}$ for each
agent, we denote by $\mathbf{a}^{j}$ the profile of actions restricted to
mechanism $j$.

  * •

Allocation rule
$\tilde{x}_{i}(\mathbf{a})=\max_{j}\tilde{x}_{i}^{j}(\mathbf{a}^{j})$. That
is, each agent is served at their highest level across all component
mechanisms.

  * •

Payment rule
$\tilde{p}_{i}(\mathbf{a})=\sum_{j}\tilde{p}_{i}^{j}(\mathbf{a}^{j})$. That
is, agents make payments to every composed mechanism.

Note that as a consequence of Definition 70, we may define the composed
feasibility environment as the set of allocation levels induced by the
component mechanisms, i.e.

$\mathcal{X}=\\{(\max_{j}x_{1}^{j},\ldots,\max_{j}x_{n}^{j})\,|\,\mathbf{x}^{1},\ldots,\mathbf{x}^{m}\in\mathcal{X}^{1},\ldots,\mathcal{X}^{m}\\}.$

The main result of this section is that weak competitive efficiency of the
component mechanisms implies weak competitive efficiency for their
simultaneous composition. That is, if each component mechanism efficiently
translates competition, measured via threshold costs, into revenue, then the
simultaneous composition of those mechanisms will do so as well.

Competitive efficiency holds with respect to cost functions $\beta_{i}$ for
the composed mechanism. for each $i$, which we may extend to the component
mechanisms in the following way. Let $a_{i}^{j}$ be an action for agent $i$ in
mechanism $j$, and let $\hat{a}_{i}^{j}$ denote the action for agent $i$ in
the composed mechanism where $i$ plays $a_{i}^{j}$ in mechanism $M^{j}$ and
withdraws from all other mechanisms $j^{\prime}\neq j$. Let $a_{i}$ be an
action in the composed mechanism consisting of actions $a_{i}^{j}$ in
mechanism $M^{j}$ for each $j$. Define
$\beta_{i}^{j}(a_{i})=\beta_{i}(\hat{a}_{i}^{j})$, where we omit conditioning
on $v_{i}$ due to the independence of $\mathcal{F}$. For any allocation
$x\in[0,1]$, we may also define agent $i$’s threshold cost
$\tau_{i}^{j}(x)=\inf_{a_{i}:\tilde{x}_{i}(a_{i})\geq x}\beta_{i}^{j}(a_{i})$,
and similarly define the threshold cost with discounted reserve
$\tau_{i}^{j,r_{i}}(x)$ and generalized threshold surplus with discounted
reserve ${\mathscr{T}}_{i}^{j,r_{i}}(x)$ for mechanism $M^{j}$ analogously to
Definitions 58 and 59. For the example of price per unit costs
$\beta_{i}^{\text{PPU}}$, $\beta_{i}^{j}$ corresponds to the price per unit
cost in the component mechanism $M^{j}$ as one would expect. We may now state
the section’s main result.

######  Theorem 71.

Let $M$ be a simultaneous composition of mechanisms $M^{1},\ldots M^{m}$. If
$M^{1},\ldots,M^{m}$ all have weak competitive efficiency $\mu$ with minimum
costs $\mathbf{r}$, then so does $M$.

The theorem will follow from two observations. First, agents’ threshold prices
are lower in the composition of mechanisms than in any individual component
mechanism. In other words, it is easier for agents to secure allocation with
more mechanisms to participate in. Second, the aggregate revenue is the sum of
the component mechanisms’ revenues. The following lemma formalizes the first
observation.

######  Lemma 72.

For any distribution of actions and values $\mathcal{G}$, any $z\in[0,1]$, any
component mechanism $M^{j}$, and any minimum cost $r_{i}$,
${\mathscr{T}}_{i}^{r_{i}}(z)\leq{\mathscr{T}}_{i}^{j,r_{i}}(z)$.

###### Proof.

By definition, $\tau_{i}(x)=\inf_{a_{i}:\tilde{x}_{i}(a_{i})\geq
x}\beta(a_{i})\leq\inf_{a_{i}:\tilde{x}_{i}(a_{i})\geq
x}\beta_{i}^{j}(a_{i})=\tau_{i}^{j}(x)$, and so
$\tau_{i}^{r_{i}}(x)\leq\tau_{i}^{j,r_{i}}(x)$. Integrating both quantities
yields ${\mathscr{T}}_{i}^{r_{i}}(x)\leq{\mathscr{T}}_{i}^{j,r_{i}}(x)$. ∎

###### Proof of Theorem 71.

Let $\mathcal{B}$ be an action distribution and $\mathcal{F}$ a product prior
over values, with joint distribution $\mathcal{G}$. Further let $\mathbf{y}$
be a feasible allocation in the composed environment, where
$\mathbf{y}^{1}\,\ldots,\mathbf{y}^{m}\in\mathcal{X}^{1},\ldots,\mathcal{X}^{m}$
denotes a profile of allocations for each component mechanism certifying the
feasibility of $\mathbf{y}$, i.e. $y_{i}=\max_{j}y_{i}^{j}$ for all $i$. For
each agent $i$, let $\delta_{ij}$ be an indicator taking value $1$ if $j$ is
the lowest index such that $y_{i}=y_{i}^{j}$, and $0$ otherwise. Further let
$\textsc{Rev}^{j}(M,\mathcal{G})$ denote the revenue from component mechanism
$M^{j}$ under $\mathcal{G}$. We obtain the following sequence of inequalities,
explained after their statement:

$\displaystyle\tfrac{1}{\mu}\textsc{Rev}(M,\mathcal{G})$
$\displaystyle=\sum\nolimits_{j}\tfrac{1}{\mu}\textsc{Rev}^{j}(M,\mathcal{G})$
$\displaystyle\geq\sum\nolimits_{j}\sum\nolimits_{i}{\mathscr{T}}_{i}^{j,r_{i}}(y_{i}^{j})$
$\displaystyle\geq\sum\nolimits_{j}\sum\nolimits_{i}{\mathscr{T}}_{i}^{r_{i}}(y_{i}^{j})$
$\displaystyle=\sum\nolimits_{j}\sum\nolimits_{i}\delta_{ij}{\mathscr{T}}_{i}^{r_{i}}(y_{i})$
$\displaystyle=\sum\nolimits_{i}{\mathscr{T}}_{i}^{r_{i}}(y_{i}).$

The first inequality follows from the definition of payments in a composed
mechanism as the sum of the revenues of the component mechanisms. The second
line comes from the assumption that each component mechanism has weak
competitive efficiency $\mu$ with minimum costs $\mathbf{r}$. The third line
follows from Lemma 72. The remaining lines follow from the definition of
feasibility in the composed mechanism. ∎

If we consider price per unit costs, then Theorem 71 implies that the
simultaneous composition of winner-pays-bid mechanisms inherits the
competitive efficiency of the component mechanisms, and hence their robust
welfare guarantees. Moreover, if all mechanisms have monopoly reserves, the
same can be said for revenue.

######  Corollary 73.

Let $M$ be the simultaneous composition of winner-pays-bid mechanisms, each
with competitive efficiency at least $\mu\leq 1$. Then in any equilibrium with
independent values satisfying no bidder communication and best response,
$\mathcal{G}=(\mathcal{F},\mathcal{B})$ with independent $\mathcal{F}$, the
welfare is at least a $\mu(e-1)/e$ fraction of optimal.

######  Corollary 74.

Let $M$ be the simultaneous composition of winner-pays-bid mechanisms, each
with competitive efficiency at least $\mu\leq 1$, and each with monopoly
reserves. Then in any equilibrium with independent values, no bidder
communication, and best response, and which respects reserves,
$\mathcal{G}=(\mathcal{F},\mathcal{B})$ with independent $\mathcal{F}$, the
revenue is at least a $\mu(e-1)/2e$ fraction of optimal.

###  5.3 All-Pay Mechanisms

This section considers single-bid all-pay mechanisms. Given an allocation rule
$\mathbf{\tilde{x}}$ mapping bids $\mathbf{b}=(b_{1},\ldots,b_{n})$ to
allocations, the corresponding all-pay mechanism has payment rule
$\tilde{p}_{i}(\mathbf{b})=b_{i}$. We will analyze the weak competitive
efficiency of all-pay mechanisms, measuring competition via gross payments
rather than price per unit of allocation. For this cost function, however,
individual efficiency fails to hold. We will instead prove that all-pay
mechanisms satisfy the weak individual efficiency of Definition 65.

Given bid distribution $\mathcal{B}$, product distribution $\mathcal{F}$ over
values, and action $a_{i}$ for agent $i$, we will consider the cost function
$\beta_{i}^{\text{AP}}(a_{i})=\tilde{p}_{i}(a_{i})=\operatorname{\mathbb{E}}_{\mathbf{a}_{-i}}\mathchoice{\left[\tilde{p}_{i}(a_{i},\mathbf{a}_{-i})\right]}{[\tilde{p}_{i}(a_{i},\mathbf{a}_{-i})]}{[\tilde{p}_{i}(a_{i},\mathbf{a}_{-i})]}{[\tilde{p}_{i}(a_{i},\mathbf{a}_{-i})]}$.
For all-pay mechanisms, payments are equal to bids, so for bid $b_{i}$, we
have $\beta_{i}^{\text{AP}}(b_{i})=b_{i}$. For such mechanisms, we therefore
also have that $\tau_{i}(x)=\inf_{b_{i}\,:\,\tilde{x}_{i}(b_{i})\geq
x}\beta_{i}^{\text{AP}}(b_{i})$ is the agent’s threshold bid for allocation
level $x$, and ${\mathscr{T}}_{i}(1)=\int_{0}^{1}\tau_{i}(x)\,dx$ is $i$’s
expected threshold bid. The connection to threshold bids will enable us to re-
use our analysis of winner-pays-bid mechanisms from Section 3. We present the
analysis in the absence of minimum costs $\mathbf{r}$, as it is typical to
consider all-pay mechanisms without reserves.

######  Theorem 75.

Let $\mathbf{\tilde{x}}$ be a deterministic allocation rule, and assume the
winner-pays-bid mechanism for $\mathbf{\tilde{x}}$ has competitive efficiency
$\mu$. Then for cost function $\beta_{i}^{\text{AP}}$, the all-pay mechanism
for $\mathbf{\tilde{x}}$ has weak competitive efficiency $\mu$.

###### Proof.

Let $\mathcal{B}$ be a bid distribution and $\mathcal{F}$ be a product
distribution over values. For a deterministic (ex post) allocation rule
$\mathbf{\tilde{x}}$, let $M^{\text{WPB}}$ and $M^{\text{AP}}$ be mechanisms
obtained by pairing $\mathbf{\tilde{x}}$ with winner-pays-bid and all-pay
payment formats, respectively. Assume $M^{\text{WPB}}$ has competitive
efficiency $\mu$. The result will follow from noting that the all-pay revenue
under $\mathcal{B}$ and $\mathcal{F}$ is weakly larger, while the threshold
surplus in $M^{\text{WPB}}$ equals the generalized threshold surplus in
$M^{\text{AP}}$. Note that we need only consider $\mathbf{y}\in\\{0,1\\}^{n}$
in the threshold surplus, as the environment is deterministic.

The revenue in a winner-pays-bid mechanism is the bids of the winners. The
revenue in an all-pay mechanism is the bids of all agents. Therefore for the
same bid distribution and prior $\mathcal{B}$ and $\mathcal{F}$,
$\textsc{Rev}(M^{\text{AP}},\mathcal{B},\mathcal{F})\geq\textsc{Rev}(M^{\text{WPB}},\mathcal{B},\mathcal{F})$.
To compare threshold surplus, note that the winner-pays-bid threshold surplus
of $(\mathcal{B},\mathcal{F})$ in $M^{\text{WPB}}$ is
$T_{i}(y_{i})=\int_{0}^{y_{i}}t_{i}(x)\,dx$. For $y_{i}=1$, this integral is
the expected threshold bid for agent $i$. Now consider the same
$(\mathcal{B},\mathcal{F})$ in $M^{\text{AP}}$. The generalized threshold
surplus for $\beta_{i}^{\text{AP}}$ is
${\mathscr{T}}_{i}(y_{i})=\int_{0}^{y_{i}}\tau_{i}(x)\,dx$. This is again
equal to $i$’s expected threshold bid in the all-pay mechanism for $y_{i}=1$.
Hence, $T_{i}(y_{i})={\mathscr{T}}_{i}(y_{i})$. We therefore have:

$\textsc{Rev}(M^{\text{AP}},\mathcal{B},\mathcal{F})\geq\textsc{Rev}(M^{\text{WPB}},\mathcal{B},\mathcal{F})\geq\mu\sum
T_{i}(y_{i})=\mu\sum{\mathscr{T}}_{i}(y_{i}),$

where the second inequality follows from the competitive efficiency of
$M^{\text{WPB}}$. ∎

We now derive a weak individual efficiency guarantee:

######  Lemma 76 (Weak Individual Efficiency of All-Pay Mechanisms).

Let $(v,b^{\theta},\tilde{x}^{\theta},\tilde{p}^{\theta})$ be a single-agent
action outcome with all-pay payments: $\tilde{p}^{\theta}(b)=b$ for all
$\theta$. If $(v,b^{\theta},\tilde{x}^{\theta},\tilde{p}^{\theta})$ satisfies
best response, then with cost function $\beta_{i}^{\text{AP}}$ the following
holds:

$\mathbb{E}_{\theta}[v\tilde{x}^{\theta}(b^{\theta})-b^{\theta}]+{\mathscr{T}}(1)\geq
v/2.$ (14)

To see that all-pay mechanisms do not satisfy the stronger definition of
individual efficiency, consider the single-bid rule $\tilde{x}(b)=b$ and an
agent with value $v=1$. Since $\tilde{u}(b)=0$ for all best responses $b$, it
follows that $\tilde{u}(b)+{\mathscr{T}}(z)\rightarrow 0$ as $z^{2}/2$, which
implies that no constant individual efficiency guarantee holds for arbitrary
$z$.

###### Proof of Lemma 76.

Let
$\tilde{u}(b^{\theta},\tilde{x}^{\theta},\tilde{p}^{\theta})=(v,\mathcal{D})$
be given. For any allocation probability $x\in[0,1]$, agent $i$ could choose
to get allocation probability at least $x$ and pay at most $\tau(x)$. Hence,
$u(\mathcal{D})=E_{\theta}[v\tilde{x}^{\theta}(b^{\theta})-b^{\theta}]\geq
vx-\tau(x)$ for all $x\in[0,1]$. We may rearrange this as $\tau(x)\geq
vx-u(\mathcal{D})$. Since we also have $t(x)\geq 0$, we may write:

$\displaystyle{\mathscr{T}}(1)\geq\int_{0}^{1}\max(vx-u(\mathcal{D}),0)\,dx$
$\displaystyle=\int_{u(\mathcal{D})/v}^{1}vx-u(\mathcal{D})\,dx$
$\displaystyle=\tfrac{v}{2}-u(\mathcal{D})+\tfrac{u(\mathcal{D})^{2}}{2v}.$

We therefore have $u(\mathcal{D})+{\mathscr{T}}(1)\geq
v/2+u(\mathcal{D})^{2}/2v$. Holding $v$ fixed and minimizing the righthand
side as a function of $u(\mathcal{D})$ yields a lower bound of $v/2$, as
desired. ∎

From Theorem 68, we may conclude that all winner-pays-bid mechanisms for
deterministic environments which have competitive efficiency $\mu$ have an
analogous all-pay mechanism with the same weak competitive efficiency. The
welfare in these mechanisms is consequently a $\mu/2$-fraction of optimal.
Notably, this includes the simultaneous composition of all-pay mechanisms.

######  Corollary 77.

For a deterministic allocation rule $\mathbf{\tilde{x}}$, let $M^{\text{WPB}}$
and $M^{\text{AP}}$ be the respective winner-pays-bid and all-pay mechanisms
for $\mathbf{\tilde{x}}$. Then if $M^{\text{WPB}}$ has competitive efficiency
$\mu\leq 1$, then in any equilibrium $(\mathcal{B},\mathcal{F})$ with product
distribution $\mathcal{F}$ satisfying no bidder communication and best
response in $M^{\text{AP}}$, the welfare is at least a $\mu/2$ fraction of
optimal.

## References

  * Allouah and Besbes, (2020) Allouah, A. and Besbes, O. (2020).  Prior-independent optimal auctions.  Management Science, 66(10):4417–4432.
  * Bergemann et al., (2017) Bergemann, D., Brooks, B., and Morris, S. (2017).  First-price auctions with general information structures: Implications for bidding and revenue.  Econometrica, 85(1):107–143.
  * Bergemann et al., (2019) Bergemann, D., Brooks, B., and Morris, S. (2019).  Revenue guarantee equivalence.  American Economic Review, 109(5):1911–29.
  * Bergemann and Morris, (2016) Bergemann, D. and Morris, S. (2016).  Bayes correlated equilibrium and the comparison of information structures in games.  Theoretical Economics, 11(2):487–522.
  * Brooks and Du, (2021) Brooks, B. and Du, S. (2021).  Optimal auction design with common values: An informationally robust approach.  Econometrica, 89(3):1313–1360.
  * Carroll, (2017) Carroll, G. (2017).  Robustness and separation in multidimensional screening.  Econometrica, 85(2):453–488.
  * Devanur et al., (2015) Devanur, N. R., Hartline, J. D., and Yan, Q. (2015).  Envy freedom and prior-free mechanism design.  Journal of Economic Theory, 156:103–143.
  * Dhangwatnotai et al., (2010) Dhangwatnotai, P., Roughgarden, T., and Yan, Q. (2010).  Revenue maximization with a single sample.  In ACM Conference on Electronic Commerce, pages 129–138.
  * Dhangwatnotai et al., (2015) Dhangwatnotai, P., Roughgarden, T., and Yan, Q. (2015).  Revenue maximization with a single sample.  Games and Economic Behavior, 91:318–333.
  * Dütting and Kesselheim, (2015) Dütting, P. and Kesselheim, T. (2015).  Algorithms against anarchy: Understanding non-truthful mechanisms.  In 16th ACM Conference on Economics and Computation.
  * (11) Edelman, B., Ostrovsky, M., and Schwarz, M. (2007a).  Internet advertising and the generalized second-price auction: Selling billions of dollars worth of keywords.  American economic review, 97(1):242–259.
  * (12) Edelman, B., Ostrovsky, M., and Schwarz, M. (2007b).  Internet advertising and the generalized second-price auction: Selling billions of dollars worth of keywords.  American economic review, 97(1):242–259.
  * Feldman et al., (2013) Feldman, M., Fu, H., Gravin, N., and Lucier, B. (2013).  Simultaneous auctions are (almost) efficient.  In Proceedings of the forty-fifth annual ACM symposium on Theory of computing, pages 201–210.
  * Forges, (1986) Forges, F. (1986).  An approach to communication equilibria.  Econometrica: Journal of the Econometric Society, pages 1375–1385.
  * Fu et al., (2015) Fu, H., Immorlica, N., Lucier, B., and Strack, P. (2015).  Randomization beats second price as a prior-independent auction.  In Proceedings of the Sixteenth ACM Conference on Economics and Computation, pages 323–323.
  * Hartline et al., (2014) Hartline, J., Hoy, D., and Taggart, S. (2014).  Price of anarchy for auction revenue.  In Proceedings of the fifteenth ACM conference on Economics and computation, pages 693–710.
  * Hartline et al., (2020) Hartline, J., Johnsen, A., and Li, Y. (2020).  Benchmark design and prior-independent optimization.  In 2020 IEEE 61st Annual Symposium on Foundations of Computer Science (FOCS), pages 294–305. IEEE.
  * Hartline et al., (2015) Hartline, J., Syrgkanis, V., and Tardos, E. (2015).  No-regret learning in bayesian games.  Advances in Neural Information Processing Systems, 28.
  * Hartline, (2013) Hartline, J. D. (2013).  Mechanism design and approximation.  Book draft. October, 122:1.
  * Hartline and Roughgarden, (2009) Hartline, J. D. and Roughgarden, T. (2009).  Simple versus optimal mechanisms.  In ACM Conference on Electronic Commerce, pages 225–234.
  * Hoy et al., (2015) Hoy, D., Nekipelov, D., and Syrgkanis, V. (2015).  Robust data-driven efficiency guarantees in auctions.  In EC Workshop on Algorithmic Game Theory and Data Science.
  * Jin and Lu, (2022) Jin, Y. and Lu, P. (2022).  First price auction is 1–1/e2 efficient.  In 2022 IEEE 63rd Annual Symposium on Foundations of Computer Science (FOCS), pages 179–187. IEEE Computer Society.
  * Kaplan and Zamir, (2012) Kaplan, T. R. and Zamir, S. (2012).  Asymmetric first-price auctions with uniform distributions: analytic solutions to the general case.  Economic Theory, 50(2):269–302.
  * Lehmann et al., (2002) Lehmann, D., Oćallaghan, L. I., and Shoham, Y. (2002).  Truth revelation in approximately efficient combinatorial auctions.  Journal of the ACM (JACM), 49(5):577–602.
  * Myerson, (1981) Myerson, R. (1981).  Optimal auction design.  Mathematics of Operations Research, 6(1):58–73.
  * Roughgarden, (2009) Roughgarden, T. (2009).  Intrinsic robustness of the price of anarchy.  In ACM Symposium on Theory of Computing, pages 513–522.
  * Roughgarden et al., (2012) Roughgarden, T., Talgam-Cohen, I., and Yan, Q. (2012).  Supply-limiting mechanisms.  In ACM Conference on Electronic Commerce, pages 844–861.
  * Syrgkanis, (2014) Syrgkanis, V. (2014).  Efficiency of mechanisms in complex markets.  PhD thesis, Cornell University.
  * Syrgkanis and Tardos, (2013) Syrgkanis, V. and Tardos, E. (2013).  Composable and efficient mechanisms.  In ACM Symposium on Theory of Computing, pages 211–220.
  * Varian, (2007) Varian, H. R. (2007).  Position auctions.  international Journal of industrial Organization, 25(6):1163–1178.