We introduce a new example of unital commutative $n$-dimensional group algebra $\R_n$ for $n \geq 2$.
The algebra $\R_n$ and the complex numbers $\C$ are astonishingly alike.
The zero divisor set of the algebra has Lebesgue $\mu_n$-measure zero.
The formula for the Haar measure is established.
Also, the analytic function theory in $\R_n,$ for $n=2k$ that similar to the classical theory in $\C$
is introduced.
This includes the Cauchy-Riemann equations, mean-value theorem and Louisville theorem.

[2000]MSC12E99 MSC13A99

§ INTRODUCTION

In this paper, we introduce a new example of unital commutative $n$-dimensional group algebra $\R_n$ for $n \geq 2$.
When $n=2$, the result is the complex numbers system $\C$.
The algebra $\R_n$ and the complex numbers $\C$ are astonishingly alike.
We name $\R_n$ the Euclidean algebra because the definition of new multiplication
works for all Euclidean space over $\R$ in a natural way.
The algebra has zero divisors  of Lebesgue $\mu_n$-measure zero.
The formula for the Haar measure is established.
We introduced the analytic function theory in $\R_n,$ for $n=2k$ that similar to the classical theory in $\C$.
This includes the Cauchy-Riemann equations, mean-value theorem and Louisville theorem.
In Section <ref> We compare $\R_n$ with other commutative hypercomplex number systems in literatures.

§ MULTIPLICATION IN $\R^N$

Let $n\in\N, n\geq 2$ and $\R^n$ be the $n$-dimensional Euclidean space over $\R$ and $\{\ve_1,\cdots, \ve_n\}$ be the standard basis of $\R^n$.
Define a mapping $\mathfrak{g}:\{\ve_1,\cdots, \ve_n\}\rightarrow \R^n$ as

\begin{align*}
\left\{
\begin{array}{llc}
\mathfrak{g}\ve_1 & =&-\ve_{n}\\
\mathfrak{g}\ve_2 & = &\ve_1\\
\cdots & \cdots & \\
\mathfrak{g}\ve_{n} &= &\ve_{n-1}\\
\end{array}
\right.
\end{align*}

Let $I_n$ denote the identity mapping (matrix) on the real $n$-dimensional Euclidean space $\R^n.$
It is clear that $\mathfrak{g}^n=-I_n$ and $\mathfrak{g}^{2n}=I_n$.
The set $\{\mathfrak{g}^\ell \mid \ell =1,2,\cdots, 2n\}$ is cyclic, hence it is an abelian group of order $2n$.
We write $\mathfrak{g}^0 = \mathfrak{g}^{2n} = I_n$.
For $0<\ell\leq n-1$, the matrix form for the mapping $\mathfrak{g}^\ell$ is

\begin{align*}
\left[\begin{array}{cc}
O & I_{n-\ell}  \\
-I_\ell & O'
\end{array}\right]
\end{align*}

where $I_{n-\ell}$ is the $(n-\ell)\times (n-\ell)$ identity matrix, $I_\ell$ is the $\ell\times \ell$ identity matrix
and $O$ a zero matrix.

Since $\mathfrak{g}^{n+k}= - \mathfrak{g}^{k}$.
The linear span of the group $\mathfrak{G}_n\equiv\{\mathfrak{g}^\ell \mid \ell =1,2,\cdots, 2n\}$ over $\R$ has dimension $n$.
Let $\M_n$ denote this linear span.
The set $\{\mathfrak{g}^0,\mathfrak{g}^1,\cdots, \mathfrak{g}^{n-1}\}$ is a basis for $\M_n$,
a linear subspace of the $n\times n$ real matrix space.
It is clear that $\M_n$ is closed under matrix multiplication.
Since the group $\{\mathfrak{g}^\ell \mid \ell =1,2,\cdots, 2n\}$ is abelian, the space $\M_n$ with the matrix multiplication is a commutative algebra of dimension $n$
with multiplicative identity $I_n$. Let $A$ be an element in $\M_n$ with $\det A \neq 0$. The inverse matrix $A^{-1}$ exists.
By the Cayley-Hamilton theorem, $A^{-1}$
is a linear combination of $\{I_n, A,A^2,\cdots,A^{n-1}\}$, hence $A^{-1}$ is also an element of $\M_n$.

Define the mapping $\varsigma:$

\begin{align*}
\varsigma : \ve_\ell \mapsto \mathfrak{g}^{\ell-1}, \ell=1,2,\cdots, n.
\end{align*}

The mapping $\varsigma$ is a bijection between bases of $\R^n$ and $\M_n$.
The mapping $\varsigma$ extends to an isomorphism from $\R^n$ onto $\M_n$ as linear spaces over $\R$.
For an element $\vu =[u_1,u_2,\cdots, u_{n-1},u_n] \in \R^n$, we have

\begin{align*}
\varsigma(\vu) = \varsigma\big(\sum_{\ell=1}^{n} u_\ell \ve_\ell\big) = \sum_{\ell=1}^{n} u_\ell \mathfrak{g}^{\ell-1}\in \M_n.
\end{align*}

For $\vu,\vv \in \R^n$, define $\vu\circledast\vv$ as

\begin{align*}
\vu\circledast\vv \equiv \varsigma^{-1}\big(\varsigma (\vu) \varsigma(\vv)\big).
\end{align*}

The Euclidean space $\R^n, n \geq 2$ equipped with the multiplication $\circledast$ is a commutative algebra with multiplicative identity $\ve_1$.

Let $\vu,\vv,\vw$ be elements in $\R^n$.

1. Associative Property
By definition we have $\varsigma(\vu\circledast\vv) = \varsigma (\vu) \varsigma(\vv)$ and
$\varsigma(\vv\circledast\vw) = \varsigma (\vv) \varsigma(\vw)$. Also the multiplication in $\M_n$ satisfies the
Associative Law. We have

\begin{align*}
(\vu\circledast\vv) \circledast \vw &= \varsigma^{-1}\big(\varsigma (\vu\circledast \vv) \varsigma(\vw)\big)\\
&= \varsigma^{-1}\big((\varsigma (\vu)\varsigma(\vv)) \varsigma(\vw)\big)\\
&= \varsigma^{-1}\big(\varsigma (\vu)(\varsigma(\vv) \varsigma(\vw))\big)\\
&= \varsigma^{-1}\big(\varsigma (\vu) \varsigma(\vv \circledast \vw)\big)\\
&=     \vu\circledast(\vv \circledast \vw).
\end{align*}

\begin{align*}
(\vu\circledast\vv) \circledast \vw = \vu\circledast(\vv \circledast \vw).
\end{align*}

2. Commutative Property
Since the multiplication in $\M_n$ is commutative, $\varsigma (\vu) \varsigma(\vv)=\varsigma (\vv) \varsigma(\vu)$, we have

\begin{align*}
\vu\circledast\vv = \varsigma^{-1}\big(\varsigma (\vu) \varsigma(\vv)\big) =
\varsigma^{-1}\big(\varsigma (\vv) \varsigma(\vu)\big) = \vv\circledast\vu .
\end{align*}

3. Distributive Property
Based on the facts that $\varsigma$ and $\varsigma^{-1}$ are isomorphisms of $\R^n$ and $\M_n$ as linear spaces over $\R$,
also the fact that $\M_n$ is an algebra satisfying the Distributive Law, we have

\begin{align*}
\vw \circledast (\vu+\vv) &= \varsigma^{-1}\big(\varsigma (\vw) \varsigma(\vu+\vv)\big)\\
&= \varsigma^{-1}\big(\varsigma (\vw) \big(\varsigma(\vu)+\varsigma(\vv)\big)\big)\\
&= \varsigma^{-1}\big(\varsigma (\vw) \varsigma(\vu)+\varsigma (\vw)  \varsigma(\vv)\big)\\
&= \varsigma^{-1}\big(\varsigma (\vw) \varsigma(\vu)\big)+\varsigma^{-1}\big(\varsigma (\vw) \varsigma(\vv)\big)\\
&= \vw \circledast \vu+\vw \circledast \vv.
\end{align*}

\begin{align*}
\vw \circledast (\vu+\vv) = \vw \circledast \vu+\vw \circledast \vv.
\end{align*}

4. Multiplicative Identity
The multiplicative identity in $\M_n$ is $I_n=\mathfrak{g}^0$ and the corresponding element in $\R^n$ is $\varsigma^{-1}(\mathfrak{g}^0)=\ve_1=[1,0,\cdots,0]$.
We have

\begin{align*}
\vu \circledast \ve_1 = \ve_1 \circledast \vu = \varsigma^{-1} \big( I_n  \varsigma (\vu) \big)=\varsigma^{-1} \big( \varsigma (\vu) \big)=\vu.
\end{align*}

The Theorem <ref> is proved.

We call this algebra the Real Euclidian Algebra. It is denoted by $\R_n$.
The multiplication table for $\circledast$ is

\begin{equation}\label{eqn:mt}
\begin{tabular}{|l||l|l|l|l|l|l|l|}
\hline
% after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
$\circledast$ & $\ve_1$ & $\ve_2$ & $\ve_3$ &... & $\ve_{n-2}$& $\ve_{n-1}$ & $\ve_n$ \\
\hline
\hline
$\ve_1$ & $\ve_1$ &$\ve_2$& $\ve_3$ &... & $\ve_{n-2}$& $\ve_{n-1}$& $\ve_n$ \\
$\ve_2$ & $\ve_2$ &$\ve_3$ & $\ve_4$ &...& $\ve_{n-1}$& $\ve_n$& -$\ve_1$ \\
$\ve_3$ & $\ve_3$ &$\ve_4$ & $\ve_5$ &...& $\ve_n$& -$\ve_1$& -$\ve_2$ \\
... & ...& ...& ...& ... &  ... &  ...&  ...\\
$\ve_{n-2}$ & $\ve_{n-2}$  & $\ve_{n-1}$ & $\ve_n$ &...& -$\ve_{n-5}$& -$\ve_{n-4}$& -$\ve_{n-3}$  \\
$\ve_{n-1}$ & $\ve_{n-1}$  & $\ve_n$ & -$\ve_1$ &...& -$\ve_{n-4}$& -$\ve_{n-3}$& -$\ve_{n-2}$  \\
$\ve_n$ & $\ve_n$ & -$\ve_1$ &-$\ve_2$ & ...& -$\ve_{n-3}$& -$\ve_{n-2}$ & -$\ve_{n-1}$ \\
\hline
\end{tabular}
\end{equation}

For $\vx,\vy\in\R_n$, if there is no confusion in the context,
we write $\vx\vy$ for $\vx\circledast\vy$ and $\vx^m$ for the product of $m$-multipliers $\vx\circledast \vx \circledast \cdots \circledast \vx$.
We reserve the notation $\vx \cdot \vy$ for the dot product. For $\va_m, \cdots \va_1,\va_0 \in \R_n$ the polynomial is well defined:

\begin{align*}
p(\vx)=\va_{m} \vx^{m} +\va_{m-1} \vx^{m-1} +\cdots + \va_{1} \vx +\va_{0}.
\end{align*}

§ REMARKS

We have some notes in order.

* The algebra $\R_n$ is the Euclidean space $\R^n$ with the standard basis $\{\ve_k \mid k=1,\cdots n\}$ and the basis follows the multiplication table (<ref>).
In particular, $\ve_1$ is the multiplicative identity $1$.
As a Euclidean space over $\R$,
for $\vx =[x_1,\cdots, x_n]$ and $\va =[a_1,\cdots, a_n]$ in $\R_n$, the distance is

\begin{align*}
|\vx-\va|=\sqrt{\sum_{j=1}^n (x_j-a_j)^2}.
\end{align*}

Also, the Lebesgue measure $\mu_n$ is naturally defined in $\R_n$.

* Let $\vx\in\R_n$, define the norm of $\|\vx\|$ as

\begin{align*}
\|\vx\| \equiv \big| \det \varsigma(\vx)  \big|.
\end{align*}

Then by the matrix multiplication properties, for $\vx,\vy\in\R_n$

\begin{align*}
\|\vx \circledast\vy\| = \big| \det \varsigma (\vx)\det \varsigma(\vy) \big|= \|\vx\| \|\vy\|
\end{align*}

* For $n\geq 2$, the group $\mathfrak{G}_n$ is

\begin{align*}
\mathfrak{G}_n=\{\mathfrak{g},\mathfrak{g}^2,\cdots, \mathfrak{g}^{n}=-1, \mathfrak{g}^{n+1}, \cdots,\mathfrak{g}^{2n}=1\}.
\end{align*}

When $n=2$,  the part between $\mathfrak{g}$ and $-1$ as well as the part between $\mathfrak{g}^{n+1}$ and $1$ reduced into empty set. So $\mathfrak{G}_2=\{\mathfrak{g},\mathfrak{g}^2,\mathfrak{g}^3,\mathfrak{g}^4\}=\{\mathfrak{g},-1,-\mathfrak{g},1\}$. This is the group
$\{i,-1,-i,1\}$ where $i$ is the imaginary unit in $\C$. So $\C$ is a natural member of $\R_n$.

* Let $\vu=[u_1,\cdots, u_n]\in \R_n$. Then

\begin{align*}
\vu&=[u_1,\cdots, u_n]\\
&=\|\vu\| \left[\frac{u_1}{\|\vu\|},\cdots, \frac{u_n}{\|\vu\|}\right]\\
&=\|\vu\| ~\big[\cos \theta_1,\cdots, \cos \theta_n \big],
\end{align*}

where $\theta_k, k=1,\cdots, n$ is the angle between $\vu$ and $\ve_k$.

For $n=2$,  we have $\varsigma(\R_2)=\M_2$.

\begin{align*}
\M_2= \left\{\left[
\begin{array}{cc}
a     & b \\
\textrm{-}b    & a \\
\end{array}
\right]
\mid a,b\in \R\right\}
\end{align*}

Also, $\det \left[
\begin{array}{cc}
a     & b \\
\textrm{-}b    & a \\
\end{array}
\right]
=a^2+b^2$. So an element $[a,b]\in \R_2$ is a zero divisor iff $a=b=0$. i.e. $\R_2$ has no non-trivial zero divisors.
Hence $\R_2$ is a $2$-dimensional commutative division algebra. By the Frobenius theorem (1878) [4]
the algebra $\R_2$ must be isomorphic to $\C$ over $\R$.

We have

\begin{align*}
\vu =[u_1,u_2] & =\sqrt{u_1^2+u_2^2} \left[\frac{u_1}{\sqrt{u_1^2+u_2^2}},\frac{u_2}{\sqrt{u_1^2+u_2^2}} \right]\\
&=\sqrt{u_1^2+u_2^2} ~\left[\cos \theta_1,\cos \theta_2 \right],
\end{align*}

Note that in the case $n=2, \theta_2=\frac{\pi}{2}-\theta_1$. This is the formula $z=|z|(\cos \theta+i\sin \theta)$ in complex numbers.

* Let $\vu,\vv\in\R_n$ and $\vv_1,\vv_2,\cdots,\vv_n$ be the $n$ row vectors in the matrix $\varsigma(\vv)^\tau$, the transpose matrix of $\varsigma(\vv)$.
It is clear we have

\begin{align*}
|\vv_j|=|\vv|, \forall j\leq n.
\end{align*}

by the multiplication table we have

\begin{align*}
\vu\vv=\vu \circledast \vv & = \big[\vu\cdot\vv_1,\cdots,\vu\cdot\vv_n\big]
= \vu\varsigma(\vv)=\vv\vu,
\end{align*}

the $j^{th}$ component of $\vu \circledast \vv$ is $\vu\cdot \vv_j$, where $``\cdot"$ is the dot product in $\R^n$.
Hence we have

Let $\vu, \vv \in\R_n$, then

\begin{align}\label{eq:linear}
\vu\circledast\vv = \vu \vv = \vv\varsigma(\vu).
\end{align}

§ ZERO DIVISORS

An element $\vx\in\R_n$ is a zero divisor if and only if $\varsigma (\vx)$ is a zero divisor in $\M_n$ if and only if
$\det(\varsigma (\vx))$,  a homogeneous polynomial in $n$ variables $\{x_1,\cdots, x_n\}$  has value zero.
Let $\mathfrak{Z}_{_{\R_n}}$ denote the set of zero divisors in $\R_n$.
For a polynomial  $p$, we will write $\mathfrak{Z}(p)$ for $\{\vx  \mid p(\vx) = 0\}$.

Let $p(\vx):\R^n\rightarrow \R$ be a non-zero polynomial on $\R^n$, then $\mu_n(\mathfrak{Z}(p)) = 0.$

Here $\mu_n$ is the Lebesgue measure on $\R^n$.

We prove the Lemma <ref> by induction on the dimensions $n$.

For $n=1$. Let $p(\vx)$ be a polynomial of degree $m\geq 1$ on $\R$. It's clear $\mathfrak{Z}(p)$ is a finite set hence has Lebesgue measure zero.

Assume the statement is true for all $n=1, \cdots, k-1$.
Consider $n=k$.

Let $p(\vx)=p(x_1,x_2,\cdots,x_k)$ be a non-zero polynomial of degree $m$ on $\R^k$ and let $t\in\R$. Define

\begin{align*}
p_{_{|x_1=t}}(\vx)\equiv p(t,x_2,\cdots,x_k).
\end{align*}

This is a function in $k-1$ variables $x_2, \cdots, x_k$.
Let $\mathcal{C}$ denote the collection of all points $t$ in $\R$ with the properties that the function $p_{_{|x_1=t}}(\vx)$  is not a zero function.
Apply the induction assumption to $p_{_{|x_1=t}}$, we have

\begin{align*}
\mu_{k-1}(\mathfrak{Z}(p_{_{|x_1=t}}))=0, \forall t \in \mathcal{C}.
\end{align*}

Let $t\in\R\backslash \mathcal{C}$.
Then $p_{_{|x_1=t}}(\vx)\equiv 0$. So $(x_1-t)$ is a factor of $p$. The polynomial has at most
finite factors. Hence the set $\R\backslash \mathcal{C}$ is a finite set. So

\begin{align*}
\mu_{k-1}(\mathfrak{Z}(p_{_{|x_1=t}}))=0, \textrm{for $t\in\R $ a.e.}
\end{align*}

By Fubini't theorem

\begin{align*}
\mu_k(\mathfrak{Z}(p))&=\int_{\R} \mu_{k-1}(\mathfrak{Z}(p_{_{|x_1=t}}))dt
= 0.
\end{align*}

The induction completes. The Lemma <ref> is proved.

Let $n\geq 3, n\in\N$ and $\mu_n$ be the Lebesgue measure on $\R^n$. Then

* The set of zero divisors in $\R_n$, $\mathfrak{Z}_{_{\R_n}}$, is a closed ideal.

* $\mu_n(\mathfrak{Z}_{_{\R_n}}) = 0$.

* The open set $\R_n\backslash \mathfrak{Z}_{_{\R_n}}$ is a group.

* $\mathfrak{Z}_{_{\R_n}}$ is star-shaped in the sense that if $\vx\in \mathfrak{Z}_{_{\R_n}}$,
then $\{t\vx \mid t \in \R\} \subset \mathfrak{Z}_{_{\R_n}}$.

Let $\G_n$ denote the group $\R_n\backslash \mathfrak{Z}_{_{\R_n}}$.

1. It is clear that $\det (\varsigma(\vu))$ is a polynomial so continuous. Therefore the set $\{\vu \mid \det(\varsigma(\vu)) \leq \frac{1}{n}\}$ is closed.
This implies that

\begin{align*}
\mathfrak{Z}_{_{\R_n}} = \bigcap_{n=1}^\infty \{\vu \mid \det(\varsigma(\vu)) \leq \frac{1}{n}\}
\end{align*}

is closed.

2. Let $\vx\in\mathfrak{Z}_{_{\R_n}}$. Then $\vx=[x_1,\cdots,x_n]$ is a zero for the polynomial $\det(\varsigma(\vx))$.
By Lemma <ref>, the conclusion follows.

3. The set $\mathfrak{Z}_{_{\R_n}}$ is closed. So its complement $\G_n$ is open. Let $\va,\vb\in\G_n$.
Then $\varsigma (\va),\varsigma (\vb) \in \M_n$, $\det \varsigma (\va) \neq 0$ and $\det \varsigma (\vb) \neq 0$.
This implies $\det (\varsigma (\va)  \varsigma (\vb)) \neq 0$, or $\va\circledast\vb\notin\mathfrak{Z}_{_{\R_n}}$.
Hence $\G_n$ is closed under multiplication $\circledast$. Notice that $\varsigma(\va)^{-1} \in \M_n$.
Let $\vc$ denote the element $\varsigma^{-1}(\varsigma(\va)^{-1})$ in $\R_n\backslash\mathfrak{Z}_{_{\R_n}}=\G_n $. We have

\begin{align*}
\va \circledast \vc = \vc \circledast \va =\varsigma^{-1}(\varsigma(\va)^{-1}  \varsigma(\va))=\varsigma^{-1}(I_n)=\ve_n.
\end{align*}

So $\vc$ is the inverse to $\va$. $\G_n$ is a group. Theorem <ref> is proved.

Zero divisors $\mathfrak{Z}_{_{\R_4}}$ in $\R_4$. Let $\vu =[u_1,u_2,u_3, u_4] \in \R^4$, then $\vu \in \mathfrak{Z}_{_{\R_4}}$
if and only if $\det(\varsigma(\vu)) = (u_1^2-u_3^2 + 2u_2u_4)^2 +(u_4^2 - u_2^2 + 2u_1u_3)^2  = 0 $ if and only if

\begin{align*}
\begin{array}{ll}
(I) \left\{
\begin{array}{l}
u_1 = s \\
u_2 = t \\
u_3 = -s + \sqrt{2}t  \\
u_4 = -\sqrt{2}s + t  \\
\end{array}
\right. & \text{ or }
(II) \left\{
\begin{array}{l}
u_1 = s \\
u_2 = t \\
u_3 = s - \sqrt{2}t \\
u_4 = \sqrt{2}s - t \\
\end{array}
\right.
\end{array}
\end{align*}

where $s$ and $t$ are free real variables.

§ THE HAAR MEASURE ON $\G_N$

Let $E$ be a Borel subset of $\G_n$ and $\mu_n$ be the Lebesgue measure on $\R^n$. Define $\nu_n (E)$ as

\begin{align*}
\nu_n(E) \equiv \int_E \frac{1}{\left|\det \varsigma(\vx)\right|} d \mu_n(\vx).
\end{align*}

$\nu_n$ is the Haar measure on $\G_n$ and it is unique up to a scalar multiplication.

Let $\va\in \G_n$ and $E$ be a Borel subset of  $\G_n$, then
$\va E = E \va \subset \G_n$.
We will show that $\nu_n(\va E)=\nu_n( E)$.
By definition,

\begin{align*}
\nu_n(\va E) &= \int_{\va E} \frac{1}{\left|\det \varsigma(\vx)\right|} d \mu_n(\vx)\\
& =\int_{E} \frac{1}{\left|\det \varsigma(\va\vy)\right|} d \mu_n(\va\vy).
\end{align*}

Lemma <ref>, we have

\begin{align*}
\va\vy&=\vy\varsigma(\va).\\
d\mu_n(\va\vy)&= |\det \varsigma(\va)| d\mu_n(\vy).
\end{align*}

\begin{align*}
\nu_n(\va E) &=\int_{E} \frac{1}{\left|\det \left(\varsigma(\va) \varsigma(\vy)\right)\right|} d \mu_n(\vy\varsigma(\va))\\
&=\int_{E} \frac{\left|\det \varsigma(\va)\right|}{\left|\det \varsigma(\va)\right| \left| \det\varsigma(\vy)\right|} d \mu_n(\vy)\\
&=\int_{E} \frac{1}{\left| \det\varsigma(\vy)\right|} d \mu_n(\vy),
\end{align*}

which is

\begin{align*}
\nu_n(\va E)=\nu_n(E).
\end{align*}

§ ANALYTIC FUNCTIONS IN $\R_{2K}$

In this section we will introduce the analytic function theory in $\R_{2k}$, which is similar to the analytic function theory in $\C$.

Let $\Omega$ be an open region in $\R_n$,
and $f$ be a function $f:\Omega\rightarrow \R_n$, that is, for a point $\vx=[x_1,\cdots,x_n]\in\Omega$,

\begin{align*}
\end{align*}

where $f_j(\vx), j=1,\cdots,n,$ are real valued functions in $n$ real variables $x_j, j=1,\cdots, n$.

Define $\vw=\ve_2 \in \R_n$ and $\vw^0 = \ve_1$. It is clear $\vw \in \G_n$ since $\det(\varsigma(\vw)) \neq 0$.
By the multiplication table, $\vw^{-j} = -\vw^{n-j}$, for $0 \leq j \leq n$.
We have

\begin{align*}
f(\vx)= \sum_{j=1} ^n f_j(\vx) \vw^{j-1}.
\end{align*}

We say a function $f(\vx)$ is differentiable in  $\Omega$ if the limit

\begin{align*}
\lim_{\vec{\delta}\rightarrow 0, \vec{\delta}\in\G_n} \frac{f(\va+\vec{\delta})-f(\va)}{\vec{\delta}}=f'(\va) \ (\textrm{ or } \frac{df}{d\vx}(\va))
\end{align*}

exists for each point $\va\in\Omega$.

In this definition, the variable  $\vec{\delta}=[x_1,x_2,\cdots,x_n]= \sum_{j=1}^n x_j\vw^{j-1}$ is approaching $0$ while remains in $\G_n$.
Notice that for $x_j \neq 0$, the vectors $x_j \vw^{j-1}$ are in $\G_n$, since $\det \varsigma(\vw^{j-1})\neq 0,  j=1,\cdots, n$.
So in the above definition of
differentiation, one can take the special path along $x_j \vw^{j-1}$ with $x_j \rightarrow 0$:

\begin{align*}
f'(\va) &= \lim_{\vec{\delta}\rightarrow 0, \vec{\delta}\in\G_n} \frac{f(\va+\vec{\delta})-f(\va)}{\vec{\delta}}
= \lim_{x_j \rightarrow 0} \frac{f(\va+x_j \vw^{j-1})-f(\va)}{x_j \vw^{j-1}},\\
f'(\va) &= \vw^{-j+1}\frac{\partial f}{\partial x_j}(\va).
\end{align*}

Let $\Omega\subset \R_n$ be an open region and $f$ be a function with domain $\Omega$.
We call $f(\vx)=[f_1(\vx),\cdots,f_n(\vx)]$ an analytic function on $\Omega$ if

* $f$ is differentiable at each point $\vx\in\Omega$.

* The second order partial derivatives of $f_m(\vx)$, $\frac{\partial^2 f_m}{\partial x_i \partial x_j}$ are continuous in $\Omega$.

It is clear $\frac{\partial^2 f_m}{\partial x_i \partial x_j}=\frac{\partial^2 f_m}{\partial x_j \partial x_i}$ by the Clairaut's Theorem.

To introduce the next results, we set $n=2k$. We have:

\begin{align}\label{eq:deri}
& f'(\va) =  \vw^{-j+1}\frac{\partial f}{\partial x_j}(\va), \text{$~j = 1, ..., 2k,$}
\end{align}

\begin{align*}
\frac{\partial f}{\partial x_{j}}(\va)
= \left[\frac{\partial f_1}{\partial x_j},\cdots,\frac{\partial f_{2k}}{\partial x_j} \right].
\end{align*}

For any given $1\leq i \leq k$, we have

\begin{align*}
&f'(\va) =  \vw^{-i+1}\frac{\partial f}{\partial x_i}(\va) = \vw^{-i-k+1}\frac{\partial f}{\partial x_{i+k}}(\va),\\
&-\vw^{k}\frac{\partial f}{\partial x_i}(\va) = \frac{\partial f}{\partial x_{i+k}}(\va),\\
O & -I_{k}  \\
I_{k} & O
\end{array}\right] \frac{\partial f}{\partial x_i}(\va)  = \frac{\partial f}{\partial x_{i+k}}(\va).
\end{align*}

Hence, if $f$ is analytic on $\Omega$, we have

\begin{align}\label{eq:laplacian}
\left\{
\begin{array}{c}
\frac{\partial f_{m}}{\partial x_i} = \frac{\partial f_{m+k}}{\partial x_{i+k}}  \\
\frac{\partial f_{m+k}}{\partial x_i} = -\frac{\partial f_{m}}{\partial x_{i+k}}~
\end{array}\right.
\forall m,i, 1\leq m,i \leq k.
\end{align}

So we have the following theorem:

Given a positive even integer $n=2k$ and an open set $\Omega\subset \R_{2k}$. For any analytic function $f=[f_1,...,f_{2k}]$ on $\Omega$, we have

* $f$ satisfies the Cauchy-Riemann equations

\begin{align*}
\left\{
\begin{array}{cc}
\frac{\partial f_m}{\partial x_i} = \frac{\partial f_{m+k}}{\partial x_{i+k}}\\
\frac{\partial f_{m+k}}{\partial x_i} = -\frac{\partial f_{m}}{\partial x_{i+k}}~
\end{array}
\right. \text{$\forall~ 1 \leq m,i \leq k$.}
\end{align*}

* For each $f_m, f_{m+k}$, $1\leq m , i\leq k$,

\begin{align}\label{eq:lap2}
\frac{\partial^2 f_m}{\partial x_i^2}+ \frac{\partial^2 f_{m}}{\partial x_{i+k}^2}= 0,
\end{align}

10pt and

\begin{align}\label{eq:lap2b}
\frac{\partial^2 f_{m+k}}{\partial x_i^2}+ \frac{\partial^2 f_{{m+k}}}{\partial x_{i+k}^2}= 0.
\end{align}

* Each $f_m$ is a harmonic function in the sense that $f_m$ has continuous second order partial derivatives and

\begin{align*}
\sum_{i=1}^{2k}\frac{\partial^2 f_m}{\partial x_i^2} = 0.
\end{align*}

We will prove the Equation (<ref>) only. The proof to the Equation (<ref>) is similar to this. Let $m \leq k$. By the Cauchy-Riemann Equations (<ref>) we have

\begin{align*}
\left\{
\begin{array}{cc}
\frac{\partial f_m}{\partial x_i} = \frac{\partial f_{m+k}}{\partial x_{i+k}}\\
~\frac{\partial f_{m}}{\partial x_{i+k}}=-\frac{\partial f_{m+k}}{\partial x_i}
\end{array}
\right. \text{$\forall~ 1 \leq i \leq k$.}
\end{align*}

\begin{align*}
\left\{
\begin{array}{cc}
\frac{\partial^2 f_m}{\partial x_i^2} = \frac{\partial^2 f_{m+k}}{\partial x_{i+k}\partial x_{i}}\\
~\frac{\partial^2 f_{m}}{\partial x_{i+k}^2}=-\frac{\partial^2 f_{m+k}}{\partial x_i\partial x_{i+k}}
\end{array}
\right. \text{$\forall~ 1 \leq i \leq k$.}
\end{align*}

\begin{align*}
\frac{\partial^2 f_m}{\partial x_i^2}+\frac{\partial^2 f_{m}}{\partial x_{i+k}^2}=
\frac{\partial^2 f_{m+k}}{\partial x_{i+k}\partial x_{i}} -  \frac{\partial^2 f_{m+k}}{\partial x_i\partial x_{i+k}}=
\frac{\partial^2 f_{m+k}}{\partial x_{i+k}\partial x_{i}} -  \frac{\partial^2 f_{m+k}}{\partial x_{i+k}\partial x_i}=0.
\end{align*}

Consider the polynomial
function in $\R_{2k}$

\begin{align*}
p(\vx)&\equiv  \va_m\vx^m +\va_{m-1}\vx^{m-1} \cdots + \va_1 \vx + \va_0\\
&=[p_1(\vx), \cdots, p_{2k}(\vx)],
\end{align*}

where $\va_m, \va_{m-1}, \cdots, \va_1,\va_0 \in\R_{2k}$.
Then each $p_j(\vx), j=1,\cdots,2k$ is a harmonic polynomial in the sense that they are polynomials and harmonic functions.

The following Example <ref> shows that, the requirement of an even dimension in Theorem <ref> is necessary.

Consider $\R_{2k-1}, k > 1$ and the function $f$ given by

\begin{align*}
f(\vx) \equiv \vx^2.
\end{align*}

The domain of $f$ is $\R_{2k-1}$.

It is left to the readers to verify that $f$ is analytic and

\begin{align*}
\end{align*}

We have

\begin{align*}
f(\vx)&=[f_1,\cdots,f_{2k-1}] = \sum_{m=1}^{2k-1} f_m \vw^{m-1} = \left(\sum_{j=1}^{2k-1} x_j\vw^{j-1}\right)^2\\
& = \sum_{i=1}^{2k-1} x_i\vw^{i-1} \sum_{j=1}^{2k-1} x_j\vw^{j-1}
= \sum_{i,j=1}^{2k-1}  x_ix_j\vw^{i+j-2}\\
f_m &=\sum_{\tiny \begin{array}{c} i+j-1=m\\ 1\leq i,j \leq 2k-1\end{array}}  x_ix_j
-\sum_{\tiny \begin{array}{c} i+j-1=m+2k-1\\1\leq i,j \leq 2k-1\end{array}}  x_ix_j.
\end{align*}

Notice that

\begin{align*}
\left(\sum_{1\leq \ell \leq 2k-1} \frac{\partial^2 }{\partial x_{\ell}^2}\right) x_ix_j
= \left\{
\begin{array}{cc}
2 & i = j;\\
0 & i \neq j.
\end{array}
\right.
\end{align*}

\begin{align*}
\lefteqn{\sum_{1\leq \ell \leq 2k-1} \frac{\partial^2 f_m }{\partial x_{\ell}^2}} \\
= &  \sum_{1\leq \ell \leq 2k-1} \frac{\partial^2 }{\partial x_{\ell}^2}
\left( \sum_{\tiny \begin{array}{c} i+j-1=m\\ 1\leq i,j \leq 2k-1\end{array}}  x_ix_j
-\sum_{\tiny \begin{array}{c} i+j-1=m+2k-1\\1\leq i,j \leq 2k-1\end{array}}  x_ix_j \right) \\
= &   \sum_{1\leq \ell \leq 2k-1} \frac{\partial^2 }{\partial x_{\ell}^2}
\left( \sum_{\tiny \begin{array}{c} 2j-1=m\\ 1\leq j \leq 2k-1\end{array}}  x_j^2
-\sum_{\tiny \begin{array}{c} 2j-1=m+2k-1\\1\leq j \leq 2k-1\end{array}}  x_j^2 \right)    \\
& \text{(Note that the set $\{j \mid 2j-1=m, 1\leq j \leq 2k-1\}$ is empty if $m$ is even.} \\
& \text{ and $\{j \mid 2j-1=m+2k-1, 1\leq j \leq 2k-1\}$ is empty if $m$ is odd. )}  \\
= & \left\{
\begin{array}{cc}
m+1 & \text{~m is odd;}\\
-m-2k & \text{~m is even.}
\end{array}
\right.
\end{align*}

Hence we have

\begin{align*}
& \sum_{\ell=1}^{2k-1}\frac{\partial^2 f_m}{\partial x_{\ell}^2} = \left\{
\begin{array}{cc}
m+1 & \text{~m is odd}\\
-m-2k & \text{~m is even}
\end{array}
\right.   \neq 0.
\end{align*}

So $f_m$ is not a harmonic function.

Next we introduce an analytic function theory in $\R_{2k}$, similar to the one in the complex number field $\C$.
We show that in two examples (Theorem <ref> and Theorem <ref>) that such a generalization is valid.

Let $f=[f_1,\cdots,f_{2k}]$ be an analytic function defined on an open region $\Omega \subset \R_{2k}$.
Let $\va\in \Omega$ and $r>0$ such that the ball $b(\va,r) \subset \Omega$. $\partial b(\va,r)$ is the boundary of the ball.
Let $\nu$ be the probability measure for the sphere.
By Theorem <ref> each $f_m, m=1,\cdots, 2k$ is a harmonic function.
It is well known that [1] for a harmonic function $u$ on $\R^n$ we have

\begin{align*}
u(\va) = \int_{\partial b(\va,r)} u(\vx) d \nu.
\end{align*}

Apply this mean value theorem to each $f_m$ we have

\begin{align*}
&= \left[\int_{\partial b(\va,r)} f_1(\vx) d \nu, \cdots, \int_{\partial b(\va,r)} f_{2k}(\vx) d \nu\right]\\
&= \int_{\partial b(\va,r)} f(\vx) d \nu.
\end{align*}

So we have Theorem <ref>, which is the Mean Value Theorem for analytic functions in $\R_{2k}$.

Let $\Omega \subset \R_{2k}$ be an open region and the ball $b(\va,r) \subset \Omega$ with center $\va\in \Omega$ and radius $r>0$.
$\partial b(\va,r)$ is the boundary of the ball and $\nu$ is the probability measure for the boundary.
For an analytic function $f=[f_1,\cdots,f_{2k}]$ defined on $\Omega$, we have

\begin{align*}
f(\va) = \int_{\partial b(\va,r)} f(\vx) d \nu.
\end{align*}

We call an analytic function an entire function on $\R_{2k}$ when it is analytic on $\Omega=\R_{2k}$.
A function $f$ is bounded in $\R_{2k}$ if there exists a positive number $M$ such that

\begin{align*}
|f(\vx)| \leq M, \forall \vx \in \R_{2k}.
\end{align*}

The next Theorem <ref> is the Louisville Theorem for Entire Functions in $\R_{2k}$.

A bounded entire function in $\R_{2k}$ is a constant function.

Let $f=[f_1, \cdots, f_{2k}]$ be an entire function on $\R_{2k}$ with a bound $M$.

\begin{align*}
|f(\vx)| \leq M, \forall \vx\in\R_{2k}.
\end{align*}

Then  for each $m\leq 2k$

\begin{align*}
|f_m (\vx)| \leq M, \forall \vx\in\R_{2k}.
\end{align*}

Notice that $f_m$ is harmonics.  By the Louisville Theorem for the harmonics functions each $f_m$ is a constant function.
Therefore $f$ is a constant function.
This proves Theorem <ref>.

§ HYPERCOMPLEX NUMBERS

In this section we will justify $\R_n$ in relation to the hypercomplex number system.
“A hypercomplex number is a number having properties departing from those of the real and complex numbers." ([9], van der Waerden (1985)).
For commutative hypercomplex number systems, the most common commutative examples are the abelian group algebra $K[G]$ where $K$ is a field and $G$ a finite abelian group, the
Davenport's $4$-dimensional commutative hypercomplex number system[3], and  the Jacobi's algebra $J_{\circledast} ^3$. [5]

The algebras $\R_3$ and $J_{\circledast} ^3$ are itentical.  Shlomo Jacobi discovered his $3$-dimensional algebra $J_{\circledast} ^3$  in 2014. His colleagues complete the unfinished manuscript after he passed away.

If we take $K=\R$ and $G=\mathfrak{G}_n$ as in Section <ref>, then the real linear span, span$(\mathfrak{G}_n)$ is the real Euclidean space $\R^n$. The group multiplication $\circledast$ will be extended
as multiplication in the linear span and the identity $\ve_1$ will be the multiplicative identity in the group algebra $\R[\mathfrak{G}_n]$ which is denoted by $\R_n$.
We name $\R_n$ the Euclidean algebra because this definition of multiplication works for all Euclidean space over $\R$ in a natural way.
In our notation, the complex number system $\C$ is $\R[\mathfrak{G}_2]=\R_2$, where the cyclic group $\mathfrak{G}_2$ is $\{i,i^2=-1,i^3=-i,i^4=1\}$.

Let $A$ be the Davenport's hypercomplex number system[3]. It has a linear basis $\{1,\vi,\vj,\vk\}$.
An element $\vu$ in $A$ is in the form

\begin{align*}
\vu = a+b\vi+c\vj+d\vk.
\end{align*}

The multiplication table is described as

\begin{align*}
&\vi \vj = \vj \vi =\vk;\\
&\vj\vk=\vk\vj=-\vi; \\
\end{align*}

In $A$, one has $\vu^2=\pm 1$ for $\vu=1,\vi,\vj,\vk$.
An element in $\R_4$ has the form

\begin{align*}
\vv = a_1\ve_1+a_2\ve_2+a_3\ve_3+a_4\ve_4.
\end{align*}

To compare $\R_4$ with $A$, we want to find all solutions that satisfies $\vv^2=\pm 1$, or equivalently:

\begin{align*}
%   \vv^2=\pm 1 ; \textrm{or equivalently }\\
(a_1\ve_1+a_2\ve_2+a_3\ve_3+a_4\ve_4)^2=\pm 1.
\end{align*}

which is

\begin{align*}
\left\{
\begin{array}{l}
a_1^2 - 2a_2a_4 - a_3^2 = \pm1\\
2a_1a_2 - 2a_3a_4= 0  \\
2a_1a_3 + a_2^2 - a_4^2 = 0 \\
2a_1a_4 + 2a_2a_3 = 0
\end{array}
\right.
\end{align*}

The solutions to the above system of equations are:

\begin{align*}
\begin{array}{ll}
\left\{
\begin{array}{l}
a_1 =\pm1 \\
a_2 = a_3 = a_4 = 0
\end{array}
\right.,
\left\{
\begin{array}{l}
a_3 =\pm1 \\
a_1 = a_2 = a_4 = 0
\end{array}
\right.,
\end{array}
\end{align*}

which means there are only $4$ solutions for $\vv$: $\{ \pm\ve_1, \pm\ve_3\}$.
However, they are linearly dependent, hence can not form a basis for $\R_4$.

It is impossible to find a new basis $\{\vv_1, \vv_1, \vv_2,\vv_4\}$ in $\R_4$ such that
$\vv_1^2 = \vv_2^2 = \vv_3^2 = \vv_4^2=\pm 1$. In other words, $\R_4$ is not algebraically isomorphic to
the Davenport's algebra $A$.

The algebra $\R_4$ is a new example of commutative hypercomplex number systerm.
We believe the following statement is true.
For each natural number $n\geq 3$, it is not possible to find a basis $\{\vb_j \mid j=1,\cdots, n\}$ for the algebra $\R_n$ with properties that

\begin{align*}
\vb_j^2=\pm1, j=1,\cdots, n.
\end{align*}

§ THE HYPERCOMPLEX NUMBER SYSTEM $\R_4$

In this last section we will list the basic facts of the $4$-dimensional hypercomplex number system $\R_4$.

* $\R_4$ is a real Euclidean space with basis $\{\ve_1,\ve_2,\ve_3, \ve_4\}$.
An element $\vv=[v_1,v_2,v_3,v_4],v_i\in\R$ in $\R_4$ has the form

\begin{align*}
\vv= v_1\ve_1 + v_2\ve_2 + v_3\ve_3 + v_4\ve_4.
\end{align*}

The distance between $\vu=[u_1,u_2,u_3,u_4]$ and $\vv=[v_1,v_2,v_3,v_4]$ is

\begin{align*}
|\vu-\vv| =  \sqrt{(u_1-v_1)^2 + (u_2-v_2)^2 + (u_3-v_3)^2+ (u_4-v_4)^2}.
\end{align*}

* The multiplication table (<ref>) for this basis becomes

\begin{equation*}
\begin{tabular}{|l||l|l|l|l|}
\hline
% after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
$\circledast$ & $\ve_1$ & $\ve_2$ & $\ve_3$ & $\ve_4$ \\
\hline
\hline
$\ve_1$ & $\ve_1$ &$\ve_2$& $\ve_3$ & $\ve_4$ \\
$\ve_2$ & $\ve_2$ &$\ve_3$ & $\ve_4$ & -$\ve_1$ \\
$\ve_3$ & $\ve_3$ &$\ve_4$ & -$\ve_1$ & -$\ve_2$ \\
$\ve_4$ & $\ve_4$ & -$\ve_1$ &-$\ve_2$ & -$\ve_3$ \\
\hline
\end{tabular}
\end{equation*}

* By Lemma <ref>

\begin{align*}
\vu\circledast\vv&=\vv\circledast\vu&=\\
&= \vv \varsigma (\vu) = \vu \varsigma(\vv)\\
&= [v_1,v_2,v_3,v_4]
\left[
\begin{array}{cccc}
u_1  &  u_2 & u_3 & u_4\\
\textrm{-}u_4 &  u_1 & u_2 & u_3\\
\textrm{-}u_3 &  \textrm{-}u_4 & u_1 & u_2\\
\textrm{-}u_2 & \textrm{-}u_3 &  \textrm{-}u_4 & u_1\\
\end{array}
\right]\\
&= [u_1,u_2,u_3,u_4]
\left[
\begin{array}{cccc}
v_1  &  v_2 & v_3 & v_4\\
\textrm{-}v_4 &  v_1 & v_2 & v_3\\
\textrm{-}v_3 &  \textrm{-}v_4 & v_1 & v_2\\
\textrm{-}v_2 & \textrm{-}v_3 &  \textrm{-}v_4 & v_1\\
\end{array}
\right]
\end{align*}

* Define

\begin{align*}
\| \vu \| = \left| \det \varsigma(\vu) \right|.
\end{align*}

\begin{align*}
\| \vu\circledast\vv\|=\| \vu\|\ \| \vv\|.
\end{align*}

* By Example <ref>, the zero divisors of $\R_4$ are the union of 2 planes:

\begin{align*}
\mathfrak{Z}_{_{\R_4}}= & \{[s,t,-s + \sqrt{2}t,-\sqrt{2}s + t] \mid s,t \in \R \} \\
& \cup \{[s,t,s - \sqrt{2}t,\sqrt{2}s - t] \mid s,t \in \R \}.
\end{align*}

* Let $E$ be a Borel subset of $\G_4=\R_4\backslash \mathfrak{Z}_{_{\R_4}}$.
By Theorem <ref>, up to a constant multiplier, the Haar measure $\nu_4$ is,

\begin{align*}
\nu_4(E) &= \int_E \frac{1}{\left|\det \varsigma(\vu)\right|} d \mu_4(\vu)\\
&=\int_E \frac{1}{(u_1^2-u_3^2 + 2u_2u_4)^2 +(u_4^2 - u_2^2 + 2u_1u_3)^2} d \mu_4,
\end{align*}

where $\mu_4$ is the Lebesgue measure in $\R^4$.

[1]
Axler, S., Bourdon, P. and Ramey, W.,
Harmonic Function Theory, Second Edition, Springer-Verlag New York, 2001

[2]
Barile, Margherita; Moslehian, Mohammad Sal; and Weisstein, Eric W.,  “Group Algebra." From MathWorld–A Wolfram Web Resource. https://mathworld.wolfram.com/GroupAlgebra.html

[3]
Davenport, C. M. A, Commutative Hypercomplex Algebra with Associated Function Theory, Clifford Algebras with Numeric and Symbolic Computations (Ed. R. Ablamowicz, P. Lounesto, and J. M. Parra). Boston, Birkh$\ddot{\textrm{u}}$ser, pp. 213-227, 1996

[4]
Frobenius, F.,
Über lineare Substitutionen und bilineare Formen, Journal f$\ddot{u}$r die reine und angewandte Mathematik 84:1-63 1878

[5]
Jacobi, S., On a novel 3D hypercomplex number system,  arXiv:1509.01459v1, 2015

[6]
Kantor, I. L.; Solodovnikov, A. S.,
Hypercomplex numbers, Springer-Verlag Berlin, New York, 1989

[7]
Olariu, S, Complex numbers in n dimensions, ELSEVIER, 2002

[8]
The Classification of Real Division Algebras, American Mathematical Monthly

[9]
van der Waerden, B. L.,  A History of Algebra from Al-Khwarizmi to Emmy Noether, New York: Springer-Verlag, 1985.


# A comparison of short-term probabilistic forecasts for the incidence of
COVID-19 using mechanistic and statistical time series models

Nicolas Banholzer  ETH Zurich, Department of Management, Technology, and
Economics, Zurich, Switzerland  University of Bern, Institute of Social and
Preventive Medicine, Bern, Switzerland  Thomas Mellan  Imperial College
London, School of Public Health, London, United Kingdom  H Juliette T Unwin
Imperial College London, School of Public Health, London, United Kingdom
Stefan Feuerriegel  LMU Munich, LMU Munich School of Management, Munich,
Germany  Munich Center for Machine Learning (MCML), Munich, Germany  Swapnil
Mishra  National University of Singapore and National University Health
System, School of Public Health and Institute of Data Science, Singapore
Samir Bhatt  Imperial College London, School of Public Health, London, United
Kingdom  University of Copenhagen, Department of Public Health, Copenhagen,
Denmark  Corresponding author<EMAIL_ADDRESS>
###### Abstract

Short-term forecasts of infectious disease spread are a critical component in
risk evaluation and public health decision making. While different models for
short-term forecasting have been developed, open questions about their
relative performance remain. Here, we compare short-term probabilistic
forecasts of popular mechanistic models based on the renewal equation with
forecasts of statistical time series models. Our empirical comparison is based
on data of the daily incidence of COVID-19 across six large US states over the
first pandemic year. We find that, on average, probabilistic forecasts from
statistical time series models are overall at least as accurate as forecasts
from mechanistic models. Moreover, statistical time series models better
capture volatility. Our findings suggest that domain knowledge, which is
integrated into mechanistic models by making assumptions about disease
dynamics, does not improve short-term forecasts of disease incidence. We note,
however, that forecasting is often only one of many objectives and thus
mechanistic models remain important, for example, to model the impact of
vaccines or the emergence of new variants.

Keywords: probabilistic forecasting, mechanistic models, time series models,
COVID-19

##  1 Introduction

Epidemic forecasting makes predictions about the development of observed
quantities such as disease incidence that can inform health policy[1, 2, 3, 4,
5, 6], e. g. by influencing decisions about resource allocation or the timing
of public health interventions to control an outbreak. Forecasts can be made
over long and short time horizons with different aims[7]. Long-term forecasts
(months to years) often analyze hypothetical scenarios or, for example, how an
epidemic develops with increased vaccination uptake or a new, more
transmissible virus strain. By contrast, short-term forecasts (days to weeks)
estimate what will happen in the near future. As a result, the performance of
short-term forecasts is usually evaluated by comparing them with the observed
outcome.

Different types of models can be used for epidemic forecasting, but these can
be very broadly categorized into mechanistic and statistical time series
models. _Mechanistic models_[8, 9] make assumptions about the data-generating
process and are based on mathematically explicit expressions about disease
dynamics. They establish a link between observed public health outcomes, such
as the reported number of cases or deaths, and unobserved outcomes directly
related to transmission, such as the number of infections or the instantaneous
reproduction number. Formalizing these links requires disease-specific
knowledge about epidemiological parameters such as the serial interval or
incubation period. _Statistical time series models_[10, 11, 12] can be
estimated without making assumptions about the data-generating process.
Instead, these models detect statistical patterns in the observed time series,
for example, through leveraging autocorrelation in the data. Since these
patterns can be detected without domain knowledge, the statistical properties
of these models are well characterized, and they can be used in a variety of
forecasting applications.

Here, we compare short-term probabilistic forecasts between mechanistic and
statistical time series models. Mechanistic models provide probabilistic
forecasts by taking into account uncertainty about the epidemiological
process. For example, mechanistic models use prior estimates of the generation
or serial interval, which are subject to uncertainty, especially during early
stages of an epidemic. Statistical time series models provide probabilistic
forecasts via Bayesian inference where model parameters are considered random
quantities to be integrated over. A comparison of the probabilistic forecasts
from mechanistic and statistical time series models can inform about the added
value of domain knowledge. If forecasts from mechanistic models were more
accurate, this would suggest that knowledge about disease dynamics and
epidemiological parameters provide added value for short-term epidemic
forecasting.

We compare the performance of short-term probabilistic forecasts from three
statistical time series models (seasonal autoregressive integrated moving
average (SARIMA)[10], Prophet[11, 13], Gaussian process (GP)[12]) with those
from two popular mechanistic models based on the renewal equation (EpiEstim[8,
14] and EpiNow2[9, 15]). We note that there are equivalences between
mechanistic models based on renewal equations and the susceptible-infected-
recovered model[16]. Our empirical comparison is based on data for the
incidence of COVID-19 from six large US states for the first pandemic year
between 15 March 2020 (the onset of the pandemic in the US) and 15 March 2021.
Models are estimated repeatedly for each day using data from the last two
months and, as in related work[7, 17, 18], probabilistic forecasts are made
for the following 2 weeks. Forecasting performance is evaluated
retrospectively by computing the difference between model-based forecasts and
observed incidence.

##  2 Materials and methods

###  2.1 Data

We use daily data on the reported number of new cases of COVID-19 for six
large US states (Arizona, California, Illinois, Maryland, New Jersey, and New
York) from 19 January 2020 to 15 March 2021. Data are retrieved using the R
package covidcast version 0.4.2[19], which aggregates data on the daily number
of confirmed cases of COVID-19 from the Johns Hopkins Coronavirus Resource
Center[20] and computes the incidence (i. e. new cases per 100,000 people)
using population data from the World Bank. Data are as of 1 January 2022, thus
including revisions to real-time data. These revisions resulted in a negative
number of new cases in a few instances, e. g. when a large number of
duplicated records were removed. We replaced negative values (which are
implausible) with the value from the day before during data preprocessing.

###  2.2 Forecasting task

The forecasting task is to make a daily 2-week forecast of the incidence of
COVID-19 at day $t$ in US state $s=1,\dots,6$. Let $N_{s,t}$ denote the
observed number of new cases of COVID-19 and let $\hat{N}_{s,t+1}$ to
$\hat{N}_{s,t+14}$ denote the daily 2-week ahead forecasts of the number of
new cases. Model-based forecasts are generated for 51 weeks (i. e. 1 year)
between 15 March 2020 and 15 March 2021. Forecasts are made on Sundays for the
following 14 days starting on Monday. Each forecast is based on a model
estimated using the historical time series up to day $t$. Specifically, the
observed data two months prior to $t$ (i. e. $N_{s,t-56},\dots,N_{s,t}$) are
used as training data for model $m$ to generate the probabilistic forecast
$\hat{N}_{s,t+k}^{m,d}$, where $\hat{N}$ is the forecasted number of new cases
and where $d=1,\dots,D$ denotes the draw (total number of $D$ draws) from the
posterior distribution $F^{m}_{s,t}$ of the model $m=1,\dots,M$ (total number
of $M$ models).

###  2.3 Summary of forecasting models

Table 1 shows a summary of the forecasting models considered in this work. The
properties of each model and the implementation are described in greater
detail in the following two sections. Note that we consider different modeling
choices for each model by varying important model parameters or
specifications. In the results, we first report for each model the modeling
choice that achieved the best overall forecasting score. The models
corresponding to these choices are then selected for comparison with each
other.

Table 1: Summary of probabilistic forecasting models and modeling choices.
Mechanistic models
---
(i) EpiEstim |  • Estimates the time-varying reproduction number over a fixed window $\tau$.
|  • Posterior distribution of reproduction number is derived analytically.
|  • Assumes cases by date of infection with a constant time shift.
|  • Uncertain generation time distribution.
|  $\Rightarrow$ Modeling choice: different values for $\tau$.
(ii) EpiNow2 |  • Extension of EpiEstim.
|  • Different methods to model the time-varying reproduction number.
|  • Uncertain delay distributions for the time from infection to reporting.
|  • Observational model with adjustment for weekday effects.
|  $\Rightarrow$ Modeling choice: different methods to model the time-varying
reproduction number.
Statistical time series models
(i) SARIMA |  • Seasonal autoregressive moving average model.
|  • Weekly differencing to account for weekday effects.
|  $\Rightarrow$ Modeling choice: different composition of the non-seasonal
part of the model.
(ii) Prophet |  • Additive time series model.
|  • Trend component with linear growth model.
|  • Growth rate changes at estimated change points determined by parameter
$\tau$.
|  • Seasonal component for weekday effects.
|  $\Rightarrow$ Modeling choice: different values for $\tau$.
(iii) GP |  • Composition of kernels to model separate aspects of the time series.
|  • Long-range squared exponential kernel for the trend.
|  • Short-range squared exponential kernel for the short-term variation.
|  • Periodic kernel for weekday effects.
|  $\Rightarrow$ Modeling choice: different length of the short-range kernel.

###  2.4 Mechanistic models

We use the following mechanistic models based on the renewal equation:
EpiEstim[8, 14] and EpiNow2[9, 15].

#### (i) EpiEstim

EpiEstim builds upon the framework by Cori et al.[8] for estimating the time-
varying reproduction number, which has been widely used for monitoring
transmission of infectious diseases[21, 22, 5, 23]. The framework assumes that
the number of new infections can be modeled with a Poisson renewal process
where the number of newly infected individuals is proportional to the number
of prevalent cases multiplied by their infectiousness[24]. Formally, the
expected incidence on day $t$ in state $s$ is estimated as

$\displaystyle E\left[I_{s,t}\right]=R_{s,t}\sum_{j=0}^{n}w_{j}I_{t-j,s},$ (1)

where $R_{s,t}$ is the time-varying reproduction number, $w$ is the generation
time (often assumed to be the serial interval distribution[25]), and
$I_{t-j,s}$ are the number of new infections in the days $j=0,\dots,n$ before
$t$. Cori et al. estimate $R_{s,t}$ with a Bayesian approach and derive the
posterior distribution of $R_{s,t}$ analytically by assuming a conjugate Gamma
prior distribution for $R_{s,t}$. The framework by Cori et al. is implemented
in the R package EpiEstim (version 2.2-4)[14].

We apply EpiEstim by assuming a constant incubation period in which case the
incidence of confirmed cases is exactly the incidence of infections, but
shifted in time. We use an uncertain generation time distribution informed by
prior estimates for the average (4.2, 95%-interval 3.3 to 5.3) and standard
deviation (4.9, 95%-interval 3.0 to 8.3) of the generation time of SARS-
CoV-2[26]. Furthermore, we assume a Gamma prior for the basic reproduction
number $R_{0}$ with a mean of 1.10 and standard deviation of 0.04 as in
related work[27]. EpiEstim estimates the reproduction number $R_{s,t}$ over a
time window $\left[t-\tau+1;t\right]$. The level of variation of the estimates
for $R_{s,t}$ can be controlled by the parameter $\tau$. Larger time windows
(i. e. higher values for $\tau$) reduce variation in the estimates of $R_{t}$
and translate into changes in incidence that are more smooth. Empirically, we
considered different modeling choices where we varied the time window from one
day to two weeks.

EpiEstim returns the estimated mean and standard deviation of $R_{t}$ in each
time window. We use the mean and standard deviation from the last time window
(i. e. the one right before the forecast is made) to compute the shape and
scale parameter for the posterior Gamma distribution of $R_{s,t}$. Based on
that, we sample posterior draws $R_{s,t}^{1},\dots,R_{s,t}^{D}$ for each
sample of the generation time distribution, which are used together with the
observed incidence $N_{s,0},\dots,N_{s,t}$ to project future incidence
$\hat{N}_{s,t+1}^{d}$ to $\hat{N}_{s,t+21}^{d}$. Incidence is modeled with a
Poisson distribution. Projections are made with the R package projections
version 0.5.4[28].

#### (ii) EpiNow2

EpiNow2 builds upon the framework by Abbott et al.[9] to estimate the time-
varying reproduction number and has been used in recent work for forecasting
the incidence of COVID-19[9, 29, 7]. EpiNow2 uses a similar approach to
EpiEstim, yet with multiple extensions. First, uncertainty about the estimated
reproduction number can be modeled using, for example, a random walk or
Gaussian process. Second, the latent number of new infections can be adjusted
to account for the proportion of the population that is still susceptible to
the virus. Third, EpiNow2 incorporates an observation model where the expected
number of new infections are mapped to the observed number of cases via
uncertain delay distributions that adjust for the time from infection to
reporting of a case. Fourth, the observed number of cases can be can be
adjusted for weekday effects. The framework by Abbott et al. is implemented in
the popular R package EpiNow2 (version 1.3.2)[15].

We apply EpiNow2 by incorporating the observational model. Uncertain delay
distributions are specified using prior estimates for the mean and standard
deviation of the generation time[26] (same as for EpiEstim), the incubation
period[30], and the reporting delay[31]. Furthermore, we assume a log-normal
prior for the basic reproduction number $R_{0}$ parameterized to have the same
mean and standard deviation as the Gamma distribution in EpiEstim.
Empirically, we considered different models to estimate the time-varying
reproduction number, using non-parametric backcalculation, a weekly random
walk, or a Gaussian process parameterized as in related work[32]. In contrast
to EpiEstim, observed cases are modeled with a negative binomial rather than a
Poisson distribution, and are adjusted for weekday effects.

###  2.5 Statistical time series models

We use the following statistical time series models: a seasonal autoregressive
integrated moving average (SARIMA) model[10], Prophet[11], and a Gaussian
process (GP)[12]. Different to mechanistic models, we directly model the log
incidence $Y_{s,t}$ (i. e. the log of the confirmed number of daily new cases
per 100,000 people). Forecasts from statistical time series models are
transformed to yield the number of number of new cases $N_{s,t}$ for
comparison with the forecasts from mechanistic models.

#### (i) SARIMA

SARIMA models[10] have been used by various studies to forecast the incidence
of COVID-19[33, 34, 35, 36, 37]. They are a linear combination of two
components: the autoregressive (AR) and the moving average (MA) component. In
the AR, the outcome $Y_{s,t}$ only depends on the lagged outcome $Y_{s,t-j}$,
where $j$ denotes the time lag. In the MA, $Y_{s,t}$ only depends on the
lagged errors $\epsilon_{s,t-j}=Y_{s,t-j}-Y_{s,t-j-1}$, i. e. the difference
between the observed and estimated outcome.

A SARIMA($p=1$,$d=0$,$q=1$)(0,1,0)7 with one AR and one MA component can be
written as

$\displaystyle\Delta Y_{s,t}=\alpha+\beta\,\Delta
Y_{s,t-1}+\phi\,\Delta\epsilon_{s,t},$ (2)

where $\Delta Y_{s,t}=Y_{s,t}-Y_{s,t-7}$ and
$\Delta\epsilon_{s,t}=\epsilon_{s,t}-\epsilon_{s,t-7}$, respectively. Note
that we use weekly differences of the log incidence in the seasonal part of
the model to account for weekday effects. Empirically, we considered different
non-seasonal parts by varying the orders of the AR ($p$), MA ($q$), and
differencing ($d$). We use weakly informative priors for the model intercept
($\alpha$), parameters of the AR ($\beta$), and the MA ($\phi$)

$\displaystyle\alpha\sim\mathrm{Normal}(0,2.5,6),\quad\beta\sim\mathrm{Normal}(0,0.5),\quad\phi\sim\mathrm{Normal}(0,0.5).$
(3)

The model is estimated with the R package bayesforecast version 1.0.1[38]

#### (ii) Prophet

Prophet[11] is a popular open source software in many forecasting applications
and has been used by various studies to model the incidence of COVID-19[39,
40, 35]. It is an additive model with three components modeling the trend of
the time series, seasonality, and holiday effects (not considered here). The
crucial component is the trend. It can be modeled assuming linear or logistic
growth and the growth rate is allowed to change at times $s_{j}$,
$j=1,\dots,J$. The changes in growth rates $\delta_{j}$ are modeled with a
hierarchical prior $\delta_{j}\sim\mathrm{Laplace}(0,\tau)$. The number of
times the growth rate can change can be controlled by adjusting the scale of
the prior ($\tau$). Greater values of $\tau$ translate into more frequent
changes of the trend.

We use Prophet assuming linear growth for the trend component and account for
weekday effects in the seasonality component. Empirically, we considered
different values for the scale of the prior for the number of change points in
the range of $\tau\in\left[0.10,0.45\right]$. The implementation is based on
the R package prophet version 1.0[13].

#### (iii) Gaussian process (GP)

GPs[12] have been used in epidemic forecasting[41, 42, 43]. A GP is a
stochastic process or collection of random variables such that the joint
distribution of every finite subset is multivariate normal (Gaussian)
distributed. A GP is completely defined by a function for the mean and
covariance (commonly called the kernel function). The mean function can be
specified arbitrarily but is often set to zero after the data is transformed,
e. g. centered. The crucial component is therefore the kernel, which
characterises the correlation between pairs of input data. Prior information
is reflected in the choice of the kernel used to describe the behavior of the
stochastic process. A popular kernel is the squared exponential (SE) kernel
$k_{\mathrm{SE}}$ defined as

$\displaystyle
k_{\mathrm{SE}}(t_{i},t_{j})=\alpha^{2}\exp\left(-\frac{\left(t_{i}-t_{j}\right)^{2}}{2\rho^{2}}\right),$
(4)

where $\alpha^{2}$ is the signal noise parameter determining the scale of the
output and $\rho$ is the length scale (or bandwidth) parameter determining the
correlation between input pairs, i. e. two time points $(t_{i},t_{j})$. Larger
$\alpha$ allow for greater deviations from the mean and larger $\rho$ result
in functions with rougher paths. Often one of the two parameters $\rho$ and
$\alpha$ is fixed due to weak identifiability[12, 44]. For time series
analysis, seasonal patterns can be modeled using a periodic kernel (PE)

$\displaystyle
k_{\mathrm{PE}}(t_{i},t_{j})=\alpha^{2}\exp\left(-\frac{2\,\sin\left(\pi\,|t_{i}-t_{j}|\,/\,p\right)^{2}}{2\rho^{2}}\right),$
(5)

where $p$ is an additional parameter for the periodicity.

We model $Y_{s,t}$ with a zero mean GP using a composition of three kernels

$\displaystyle Y_{s,t}$
$\displaystyle\sim\mathrm{GP}\left(0,k_{\mathrm{SE}}^{\mathrm{long}}+k_{\mathrm{SE}}^{\mathrm{short}}+k_{\mathrm{PE}}^{\mathrm{week}}+\sigma^{2}\delta_{ij}\right)\leavevmode\nobreak\
,$ (6)

where $k_{\mathrm{SE}}^{\mathrm{long}}$ is a SE kernel with a long length
scale and $k_{\mathrm{SE}}^{\mathrm{short}}$ is SE kernel with short length
scale, $k_{\mathrm{PE}}^{\mathrm{week}}$ is PE kernel for weekday effects, and
$\sigma^{2}\delta_{ij}$ is a regularizing term with a Kronecker delta function
ensuring $\sigma^{2}$ Gaussian noise is only added when $i=j$. The SE kernels
model the long-term trend and the short-term variation using a long- and
short-range correlation structure[44].

Prior specifications for each kernel are as follows. For
$k_{\mathrm{SE}}^{\mathrm{long}}$, we set the length scale to $\rho=56$ days
(the length of the training period) and use a weakly informative prior for
$\alpha\sim\text{Student-t}(\nu=3,\mu=0,\sigma=5)$. For
$k_{\mathrm{PE}}^{\mathrm{week}}$, we set the length scale to $\rho=1$, the
periodicity to $p=7$ days, and use
$\alpha\sim\text{Student-t}(\nu=5,\mu=0,\sigma=2)$. Empirically, we considered
different length scales $\rho\in(1,7,14)$ days in the short-range kernel
$k_{\mathrm{SE}}^{\mathrm{short}}$, allowing the GP to react more or less
quickly to short-term variation. A weakly informative prior is given by
$\alpha\sim\text{Student-t}(\nu=7,\mu=0,\sigma=2)$. We implemented the GPs in
the probabilistic programming language Stan version 2.21.0[45].

###  2.6 Evaluation measures

We evaluate probabilistic forecasts over a 2-week time window. We aggregate
daily forecasts by week to eliminate potential weekday effects, which is not
accounted for in the EpiEstim model. The following evaluation measures are
separately computed for the 1- and 2-week ahead forecast.

Forecasting performance is primarily assessed based on the continuous ranked
probability score (CRPS)[46]. It measures the difference between the
forecasted and the observed value and is a strictly proper scoring rule. It is
defined as

$\displaystyle CRPS(F_{f},x_{t})=\int\big{[}F(y_{t})-\mathbf{1}(y_{t}\geq
x_{t})\big{]}^{2}\,\mathrm{d}y_{t},$ (7)

where $F$ is the cumulative distribution function of the forecast $y$ for
observation $x$ at time $t$. Note that the CRPS is a generalization of the
absolute error between the forecasted and observed value. A lower CRPS
therefore indicates a more accurate forecast.

We follow recommendations to compute CRPS for the log of the forecasted
incidence[47]. This corresponds to evaluating the multiplicative rather than
the absolute forecasting error. Importantly, taking the log mitigates the
influence of extreme outliers in the forecasts often made by models during an
exponentially growing epidemic.

CRPS is computed for each forecast and scores are summarized by model, state,
and epidemic phase (labeled manually, see S1 Appendix Fig. ??). In addition,
we perform pairwise model comparisons[48], thereby measuring the forecasting
performance of each model relative to others. This is accomplished by first
computing the ratio of scores. The relative skill is then computed by taking
the geometric mean of all score ratios. Lower score ratios (and lower relative
skill) indicate better forecasting performance in pairwise comparisons.

Secondary performance measures include sharpness, calibration, and bias.
Sharpness is measured as the normalized median of the absolute deviation from
the median of the predictive samples[49]. Calibration is assessed based on the
density of the probability integral transform (PIT) values[46] as well as the
empirical coverage of the 50 % and 95 % predictive interval (PI)[7]. Bias is
evaluated with the average proportion of posterior draws that are above the
observed value.

Primary and secondary performance measures are computed using the R packages
scoringRules version 1.0.2[50] and scoringutils[48] version 1.0.1.

An important quality of epidemic forecasts is to anticipate exponential
growth. This quality is not captured specifically by the CRPS. Therefore, we
consider an additional metric by McDonald et al.[51] for detecting so-called
“hotspots”. A hotspot at time $t$ in state $s$ is defined as a weekly increase
in incidence of more than 25%, i. e.

$\displaystyle H_{s,t}=\begin{cases}1,\quad\text{if }\leavevmode\nobreak\
(I_{s,t}-I_{s,t-7})/(I_{s,t-7})\geq 0.25,\\\
0,\quad\text{otherwise.}\end{cases}$ (8)

The probability for a hotspot predicted by model $m$ is then equal to the
proportion of posterior draws where the weekly change in the forecast is
larger than 25%. Performance is assessed with the area under the receiving
operating characteristic curve (AUC). Similar to McDonald et al.[51], we
exclude data where the weekly incidence is below 70 new cases per 100,000
people. This removes 35% of the observations. Of the remaining observations,
25% are classified as hotspots.

###  2.7 Software and model estimation

All analyses are performed in the software R (version 4.2.0)[52].
Probabilistic models are estimated using Markov chain Monte Carlo (MCMC)
sampling (except for EpiEstim where the posterior is derived analytical). Each
model is estimated with 4 Markov chains and 1,000 iterations of which the
first 500 iterations were discarded as part of the warm-up. This results in
2,000 posterior samples.

###  2.8 Data and code availability

Data are from publicly available data sources (Johns Hopkins Coronavirus
Resource Center for epidemiological data; World Bank for population data).
Preprocessed data files together with reproducible code is available from
https://github.com/nbanho/covid_predict.

##  3 Results

The specific modeling choice has only a marginal influence on the overall
forecasting performance of each model (S1 Appendix Fig. ??). Table 2 shows the
modeling choice with the highest overall CRPS for each model. The models
corresponding to these choices are used for the subsequent between-model
comparisons of mechanistic and statistical time series models.

Table 2: Modeling choices with highest overall CRPS. Results for all choices are shown in S1 Appendix Fig. ??. Model | Modeling choice with highest overall CRPS
---|---
Mechanistic models
EpiEstim | Time-varying reproduction number estimated over a 7-day time window.
EpiNow2 | Time-varying reproduction number modeled with a Gaussian process.
Statistical time series models
SARIMA |  SARIMA(1,0,1)(0,1,0)7 model with one AR and one MA component.
Prophet |  Prior scale $\tau=0.45$ allowing for high number of change points for the growth rate.
GP |  Short-range kernel with a length scale of $\rho=7\,$days.

The 1- and 2-week ahead forecasts by each model for Arizona are shown as
examples in Fig. 1. The 1-week ahead forecast is well-calibrated across all
models and covers the observed weekly incidence most of the time, except for
the peak in January 2021 for SARIMA and Prophet. Coverage becomes worse for
the 2-week ahead forecast, especially the peak in January 2021 is
overestimated by most models and the subsequent decline is not anticipated.
Forecasts from EpiNow2 are wider than those from EpiEstim. Similarly,
forecasts from GP are wider than those from SARIMA and Prophet.

Fig 1: Probabilistic forecasts for Arizona by model. Probabilistic forecast
(posterior mean as red line, 95%-PI as shaded area) for Arizona by each model.
(a) 1-week ahead forecast. (b) 2-week ahead forecast. Observed incidence is
shown with black lines. Probabilistic forecasts for all US states by model are
shown in S1 Appendix ??-??.

Results for secondary performance measures (calibration, sharpness, and bias)
across states are shown in Fig. 2. The distribution of the probability
integral transform for EpiEstim, SARIMA, and Prophet is U-shaped (Fig. 2a),
indicating that the predictive distribution is too narrow, whereas it is hump-
shaped for EpiNow2 and GP, indicating that the predictive distribution is too
wide. GP and EpiNow2 achieve good coverage of the 50%-PI and the 95%-PI,
whereas EpiEstim, SARIMA, and Prophet cover both PIs less often than
theoretically predicted (Fig. 2b). Forecasts from EpiNow2 are least sharp,
especially the 2-week ahead forecasts exhibit by far the highest dispersion
(Fig. 2c). EpiEstim tends to overestimate and GP tends to underestimate the
true incidence (Fig. 2d). Other models do not show a similar behavior in terms
of bias.

Fig 2: Evaluation of secondary performance measures overall. (a) Probabilistic
calibration assessed with the density distribution of the probability integral
transform. (b-c) Coverage of the 50 % (left panel) and 95 % (right panel)
predictive interval [PI]. (d) Sharpness assessed with dispersion measured
median of the absolute deviation from the median of the predictive samples.
(e) Bias assessed with the average proportion of posterior draws above
observed the incidence.

CRPS scores are shown in Fig. 3. The distribution of scores are right-skewed
for all models, indicating a considerable number of large forecasting errors
(Fig. 3a). For the 1-week ahead forecast, the scores are comparable between
models. Mean and median CRPS are lowest (best) for GP and highest for Prophet.
The latter also has the highest CRPS for the 95% quantile. CRPS scores are
more dispersed for the 2-week ahead forecasts, but the differences between
models remain similar, albeit more visible. GP achieves the best scores, and
the distribution is more concentrated on low CRPS. By contrast, the
distribution is more dispersed for Prophet and the CRPS for the 90% and 95%
quantiles are higher. CRPS scores for SARIMA, EpiNow2, and EpiEstim are
comparable and in between those of GP and Prophet. Differences in forecasting
performance are confirmed by pairwise model comparisons. GP has the lowest
(best) relative skill for both the 1- and 2-week ahead forecast (Fig. 3b) .
Compared against other models, GP achieves the highest (best) score ratio,
although the difference to SARIMA and EpiEstim are statistically not
significant for the 2-week ahead forecast (Fig. 3c). In addition, GP has the
lowest CRPS scores in all six US states (S1 Appendix Fig. ??) and performs
particularly well in phases of (sub)exponential decline (S1 Appendix Fig. ??)
because it anticipates sooner the peak and subsequent decline of incidence.

Fig 3: Evaluation of forecasting performance. (a) Empirical quantile
distribution of continuous ranked probability score (CRPS) by model. Each dot
represents one of the 20 quantiles of the empirical distribution. Median CRPS
is shown with the blue square and mean CRPS with the black dot. Lower CRPS
indicate better forecasting performance. (b) Relative skill for pairwise model
comparisons according to CRPS scores. Lower values indicate better relative
forecasting performance. (c) Mean score ratio for GP compared against other
models. Score ratios closer to 1 indicate forecasting performances more
similar to GP.

The CRPS evaluates the overall forecasting performance. Another important
quality of forecasting models is to predict hotspots, i. e. large increases in
weekly incidence. A comparison of the AUC for predicting hotspots shows that
statistical time series models overall predict hotspots more accurately than
mechanistic models (Fig. 4). Prophet has the highest AUC for both the 1- and
2-week ahead forecast, followed by GP and SARIMA.

Fig 4: Evaluation of hotspot prediction. Area under the receiver operating
characteristic curve (AUC) for predicting hotspot, defined as a weekly
increase in incidence of more than 25% to more than 10 new cases per 100,000
people). Higher AUC indicate greater ability to detect hotspots.

##  4 Discussion

###  4.1 Summary

We compared the short-term forecasts from mechanistic and statistical time
series models for the incidence of COVID-19 in six large US states during the
first year of pandemic. Forecasts from statistical time series models were at
least as accurate as those from mechanistic models. GP achieved the best
overall forecasting performance, but the difference to other models such as
SARIMA and EpiEstim were not large, especially for the 2-week ahead forecast.
The forecasts from GP and EpiNow2 were less sharp than those from other models
(EpiEstim, SARIMA, and Prophet), but their empirical coverage was better.
Finally, statistical time series models predicted large increases in incidence
more accurately than mechanistic models.

###  4.2 Findings in context

Our work shows that statistical time series models achieve similar if not
better short-term forecasting performance than mechanistic models. This is in
line with previous work showing that statistical time series models with
various known applications can also provide reasonably accurate forecasts for
disease spread[35, 36, 41]. Nevertheless, epidemic forecasts remain difficult
for both types of models. The reason is that forecasts are better understood
as projections that are entirely based on the historical time series of the
outcome. As such, these projections assume that nothing is changing and that
the historical trend continues into the future. However, there can be factors
influencing disease spread that are not captured by the historical time
series, e. g. voluntary behavioral changes or non-pharmaceutical
interventions[53, 54, 25].

A related study compared human- and model-based forecasts for the cases and
deaths of COVID-19 up to four weeks ahead[32]. Human-based forecasts were
elicited using a web application where participants could select a predictive
distribution (e. g. log-normal), and then visualize and interactively adjust
their numerical forecasts. Individual forecasts were aggregated with a
quantile-wise mean and the crowd forecasts were compared to the probabilistic
forecasts from mechanistic models. The study showed that crowd forecasts were
more accurate than forecasts from mechanistic models. In particular,
mechanistic models often overestimated the incidence at the epidemic peaks.
Similarly, we found in our study that the peak incidence was often
overestimated by mechanistic models, but also statistical time series models.
An explanation for why model-based forecasts overestimate the peak compared to
human-based forecasts may be that only the latter can factor in the possible
impact of non-pharmaceutical interventions that are implemented right at the
time when forecasts have to be made.

Model-based forecasts of disease incidence could be further improved by using
ensembles or including auxiliary indicators. Previous studies showed that the
combined forecast of multiple models can be more accurate than the forecasts
from each individual model[7, 55, 56], and that additional data streams (e. g.
online surveys or mobility) can improve forecasting performance[51, 57, 58,
59]. Ensemble forecasts and auxiliary indicators were not the focus of this
work, but we acknowledge their potential.

###  4.3 Limitations

Our study has limitations. First, we evaluated the multiplicative error of
forecasts by comparing the log incidence. We chose the multiplicative error as
it is less sensitive to outliers and ensures comparability across states and
epidemic phases[32]. A drawback of the multiplicative error is that it is the
same regardless of the absolute level of incidence. This can be unfavorable to
models with large errors during phases of low incidence, which are less
concerning. Second, we used the CRPS as our primary performance measure to
evaluate forecasts. The CRPS is a proper scoring rule that is used across
multiple domains to evaluate numerical forecasts[46]. However, policymakers
may evaluate forecasts differently using a utility function which, for
example, penalizes underestimation of disease incidence more than
overestimation. Third, although we compared models that are commonly used to
forecast incidence, they still only represent a subset of all possible models.
We therefore acknowledge that superior performance may be achieved with
extensions of the models considered in this work or other types of models that
are not considered in this study, e. g. compartmental models[60], gradient
boosted trees[36], or deep learning[59]. Fourth, forecasts were made and
analyzed retrospectively, which may affect their evaluation. By contrast,
forecasts made in real-time are often based on incomplete incidence data[51]
or forecasters may make manual adjustments before submission[7]. Fifth, we
evaluated our forecasts based on the reported incidence of COVID-19. We
therefore acknowledge that the actual number of cases of COVID-19 may be
higher than reported number of cases[61], i. e. ascertainment bias is not
addressed during our analysis. Sixth, we evaluated forecasts on the incidence
of COVID-19 in six large US states over one year with carefully labeled
epidemic phases. Future research may extend our analysis, e. g. considering
different geographical regions, other epidemiological outcomes (e. g.
deaths[7, 55] or hospitalizations[37]), and or even epidemics of other
pathogens (e. g. influenza[1, 3]).

###  4.4 Concluding remarks and outlook

Mechanistic models incorporate knowledge about disease dynamics. Such domain
knowledge is important, for example, when modeling the effects of vaccination
uptake or the emergence of new virus variants. However, the added value for
short-term forecasts of disease incidence is less clear. We empirically
compared short-term probabilistic forecasts of disease from mechanistic and
statistical time series models. We found that statistical time series models
achieved similar overall performance than mechanistic models, suggesting that
disease incidence can be predicted in the short term without modeling disease
dynamics. Future research may evaluate whether forecasts can be improved
further by other means such as considering ensemble models or auxiliary
indicators.

## References

  * [1] Biggerstaff M, Johansson M, Alper D, Brooks LC, Chakraborty P, Farrow DC, et al.  Results from the second year of a collaborative effort to forecast influenza seasons in the United States.  Epidemics. 2018;24:26–33.  doi:10.1016/j.epidem.2018.02.003.
  * [2] Meltzer MI, Atkins CY, Santibanez S, Knust B, Petersen BW, Ervin ED, et al.  Estimating the future number of cases in the Ebola epidemic – Liberia and Sierra Leone, 2014–2015.  Morbidity and Mortality Weekly Report. 2014;63(3).  Available from: https://stacks.cdc.gov/view/cdc/24901.
  * [3] Bogoch II, Brady OJ, Kraemer MUG, German M, Creatore MI, Kulkarni MA, et al.  Anticipating the international spread of Zika virus from Brazil.  Lancet. 2016;387(10016):335–336.  doi:10.1016/S0140-6736(16)00080-5.
  * [4] Nsoesie E, Mararthe M, Brownstein J.  Forecasting Peaks of Seasonal Influenza Epidemics.  PLoS Curr. 2013;5.  doi:10.1371/currents.outbreaks.bb1e879a23137022ea79a8c508b030bc.
  * [5] Bhatia S, Parag KV, Wardle J, Imai N, Van Elsland SL, Lassmann B, et al.. Global predictions of short- to medium-term COVID-19 transmission trends: A retrospective assessmen.  medRxiv [Preprint]; 2021 [cited 2022 March 3].  Available from: https://www.medrxiv.org/content/10.1101/2021.07.19.21260746.abstract. doi:10.1101/2021.07.19.21260746.
  * [6] Massonnaud C, Roux J, Crépey P. COVID-19: Forecasting short term hospital needs in France.  medRxiv [Preprint]; 2020 [cited 2022 March 3].  Available from: https://www.medrxiv.org/content/10.1101/2020.03.16.20036939v1. doi:10.1101/2020.03.16.20036939.
  * [7] Bracher J, Wolffram D, Deuschel J, Görgen K, Ketterer JL, Ullrich A, et al.  A pre-registered short-term forecasting study of COVID-19 in Germany and Poland during the second wave.  Nat Commun. 2021;12:5173.  doi:10.1038/s41467-021-25207-0.
  * [8] Cori A, Ferguson NM, Fraser C, Cauchemez S.  A new framework and software to estimate time-varying reproduction numbers during epidemics.  Am J Epidemiol. 2013;178(9):1505–1512.  doi:10.1093/aje/kwt133.
  * [9] Abbott S, Hellewell J, Thompson RN, Sherratt K, Gibbs HP, Bosse NI, et al.  Estimating the time-varying reproduction number of SARS-CoV-2 using national and subnational case counts.  Welcome Open Research. 2020;(5):112.  doi:10.12688/wellcomeopenres.16006.2.
  * [10] Hyndman RJ, Athanasopoulos G.  Chapter 8: ARIMA models.  In: Forecasting: principles and practice. OTexts; 2018. p. 221–274.
  * [11] Taylor SJ, Letham B.  Forecasting at Scale.  Am Stat. 2018;72(1):37–45.  doi:10.1080/00031305.2017.1380080.
  * [12] Rasmussen CE.  Gaussian Processes for Machine Learning.  MIT Press; 2006.
  * [13] Taylor S, Letham B. Prophet: Automatic Forecasting Procedure; 2021.  R package version 1.0.  Available from: https://cran.r-project.org/web/packages/prophet.
  * [14] Cori A, Kamvar Z, Stockwin J, Jombart T, Dahlqwist E, FitzJohn R, et al.. EpiEstim: A tool to estimate time varying instantaneous reproduction number during epidemics; 2021.  R package version 2.2-4.  Available from: https://cran.r-project.org/web/packages/EpiEstim.
  * [15] Abbott S, Hellewell J, Sherrat K, Gostic KM, Hickson J, Badr HS, et al.. EpiNow2: Estimate real-time case counts and time-varying epidemiological parameters; 2021.  R package version 1.3.2.  Available from: https://cran.r-project.org/web/packages/EpiNow2.
  * [16] Champredon D, Dushoff J, Earn DJD.  Equivalence of the Erlang-distributed SEIR epidemic model and the renewal equation.  SIAM J Appl Math. 2018;78(6):3258–3278.
  * [17] Jewell NP, Lewnard JA, Jewell BL.  Predictive mathematical models of the COVID-19 pandemic: Underlying principles and value of projections.  JAMA. 2020;323(19):1893–1894.  doi:10.1001/jama.2020.6585.
  * [18] Krymova E, Béjar B, Thanou D, Sun T, Manetti E, Lee G, et al.  Trend estimation and short-term forecasting of COVID-19 cases and deaths worldwide.  Proc Natl Acad Sci U S A. 2022;119(32):e2112656119.  doi:10.1073/pnas.2112656119.
  * [19] Arnold T, Bien J, Brooks L, Colquhoun S, Farrow D, Grabman J, et al.. Covidcast: Client for delphi’s COVIDcast epidata API; 2021.  R package version 0.4.2.  Available from: https://cran.r-project.org/web/packages/covidcast.
  * [20] Dong E, Du H, Gardner L.  An interactive web-based dashboard to track COVID-19 in real time.  Lancet Infect Dis. 2020;20(5):533–534.  doi:10.1016/S1473-3099(20)30120-1.
  * [21] Riccardo F, Ajelli M, Andrianou XD, Bella A, Manso MD, Fabiani M, et al.  Epidemiological characteristics of COVID-19 cases and estimates of the reproductive numbers 1 month into the epidemic, Italy, 28 January to 31 March 2020.  Euro Surveill. 2020;25(49):2000790.  doi:10.2807/1560-7917.ES.2020.25.49.2000790.
  * [22] Al Wahaibi A, Al Manji A, Al Maani A, Al Rawahi B, Al Harthy K, Alyaquobi F, et al.  COVID-19 epidemic monitoring after non-pharmaceutical interventions: The use of time-varying reproduction number in a country with a large migrant population.  Int J Infect Dis. 2020;99:466–472.  doi:10.1016/j.ijid.2020.08.039.
  * [23] Huisman JS, Scire J, Angst DC, Li J, Neher RA, Maathuis MH, et al.  Estimation and worldwide monitoring of the effective reproductive number of SARS-CoV-2.  eLife. 2022;11:e71345.  doi:10.7554/eLife.71345.
  * [24] Fraser C.  Estimating individual and household reproduction numbers in an emerging epidemic.  PLoS One. 2007;2(8):e758.  doi:10.1371/journal.pone.0000758.
  * [25] Flaxman S, Mishra S, Gandy A, Unwin HJT, Mellan TA, Coupland H, et al.  Estimating the effects of non-pharmaceutical interventions on COVID-19 in Europe.  Nature. 2020;584(7820):257–261.  doi:10.1038/s41586-020-2405-7.
  * [26] Ganyani T, Kremer C, Chen D, Torneri A, Faes C, Wallinga J, et al.  Estimating the generation interval for coronavirus disease (COVID-19) based on symptom onset data, March 2020.  Euro Surveill. 2020;25(17):2000257.  doi:10.2807/1560-7917.ES.2020.25.17.2000257.
  * [27] Bosse NI, Abbott S, Bracher J, Hain H, Quilty BJ, Jit M, et al.. Comparing human and model-based forecasts of COVID-19 in Germany and Poland.  medRxiv [Preprint]; 2021 [cited 2022 March 3].  Available from: https://www.medrxiv.org/content/10.1101/2021.12.01.21266598v1. doi:10.1101/2021.12.01.21266598.
  * [28] Jombart T, Nouvellet P, Bhatia S, Kamvar ZN, Taylor T, Ghozzi S. Projections: Project Future Case incidence; 2021.  R package version 0.5.4.  Available from: https://cloud.r-project.org/web/packages/projections.
  * [29] Davies NG, Abbott S, Barnard RC, Jarvis CI, Kucharski AJ, Munday JD, et al.  Estimated transmissibility and impact of SARS-CoV-2 lineage B.1.1.7 in England.  Science. 2021;372(6538):eabg3055.  doi:10.1126/science.abg3055.
  * [30] Lauer SA, Grantz KH, Bi Q, Jones FK, Zheng Q, Meredith HR, et al.  The incubation period of Coronavirus Disease 2019 (COVID-19) from publicly reported confirmed cases: Estimation and application.  Ann Intern Med. 2020;172(9):577–582.  doi:10.7326/M20-0504.
  * [31] Cereda D, Tirani M, Rovida F, Demicheli V, Ajelli M, Poletti P, et al.. The early phase of the COVID-19 outbreak in Lombardy, Italy.  arXiv [Preprint]; 2020 [cited 2022 March 3].  Available from: https://arxiv.org/abs/2003.09320. doi:10.48550/ARXIV.2003.09320.
  * [32] Bosse NI, Abbott S, Bracher J, Hain H, Quilty BJ, Jit M, et al.  Comparing human and model-based forecasts of COVID-19 in Germany and Poland.  PLoS Comput Biol. 2022;18(9):e1010405.  doi:10.1371/journal.pcbi.1010405.
  * [33] Alabdulrazzaq H, Alenezi MN, Rawajfih Y, Alghannam BA, Al-Hassan AA, Al-Anzi FS.  On the accuracy of ARIMA based prediction of COVID-19 spread.  Results Phys. 2021;27:104509.  doi:10.1016/j.rinp.2021.104509.
  * [34] Kufel T.  ARIMA-based forecasting of the dynamics of confirmed COVID-19 cases for selected European countries.  Equilibrium. 2020;15(2):181–204.  doi:10.24136/eq.2020.009.
  * [35] Wang Y, Yan Z, Wang D, Yang M, Li Z, Gong X, et al.  Prediction and analysis of COVID-19 daily new cases and cumulative cases: times series forecasting and machine learning models.  BMC Infect Dis. 2022;22(1):495.  doi:10.1186/s12879-022-07472-6.
  * [36] Rahman MS, Chowdhury AH, Amrin M.  Accuracy comparison of ARIMA and XGBoost forecasting models in predicting the incidence of COVID-19 in Bangladesh.  PLoS Global Public Health. 2022;2(5):e0000495.  doi:10.1371/journal.pgph.0000495.
  * [37] Perone G.  Comparison of ARIMA, ETS, NNAR, TBATS and hybrid models to forecast the second wave of COVID-19 hospitalizations in Italy.  Eur J Health Econ. 2022;23(6):917–940.  doi:10.1007/s10198-021-01347-4.
  * [38] Matamoros AA, Torres CC, Dala A, Hyndman R, O’Hara-Wild M. Bayesforecast: Bayesian Time Series Modeling with Stan; 2021.  R package version 1.0.1.  Available from: https://cran.r-project.org/web/packages/bayesforecast.
  * [39] Battineni G, Chintalapudi N, Amenta F.  Forecasting of COVID-19 epidemic size in four high hitting nations (USA, Brazil, India and Russia) by Fb-Prophet machine learning model.  Appl Comput Inform. 2020;doi:10.1108/ACI-09-2020-0059.
  * [40] Aditya Satrio CB, Darmawan W, Nadia BU, Hanafiah N.  Time series analysis and forecasting of coronavirus disease in Indonesia using ARIMA model and PROPHET.  Procedia Comput Sci. 2021;179:524–532.  doi:10.1016/j.procs.2021.01.036.
  * [41] Johnson LR, Gramacy RB, Cohen J, Mordecai E, Murdock C, Rohr J, et al.  Phenomenological forecasting of disease incidence using heteroskedastic Gaussian processes: A dengue case study.  Ann Appl Stat. 2018;12(1):27–66.  doi:10.1214/17-AOAS1090.
  * [42] Arias Velásquez RM, Mejía Lara JV.  Forecast and evaluation of COVID-19 spreading in USA with reduced-space Gaussian process regression.  Chaos Solitons Fractals. 2020;136:109924.  doi:10.1016/j.chaos.2020.109924.
  * [43] Qian Z, Alaa AM, van der Schaar M.  In: Advances in Neural Information Processing Systems. vol. 33; 2020. p. 10729–10740.  Available from: https://proceedings.neurips.cc/paper/2020/hash/79a3308b13cd31f096d8a4a34f96b66b-Abstract.html.
  * [44] Hawryluk I, Hoeltgebaum H, Mishra S, Miscouridou X, Schnekenberg RP, Whittaker C, et al.  Gaussian process nowcasting: Application to COVID-19 mortality reporting.  Proceedings of the 37th Conference on Uncertainty in Artificial Intelligence. 2021;161:1258–1268.  Available from: https://proceedings.mlr.press/v161/hawryluk21a.html.
  * [45] Carpenter B, Gelman A, Hoffman MD, Lee D, Goodrich B, Betancourt M, et al.  Stan: A probabilistic programming language.  J Stat Softwb. 2017;76:1–32.  doi:10.18637/jss.v076.i01.
  * [46] Gneiting T, Balabdaoui F, Raftery AE.  Probabilistic forecasts, calibration and sharpness.  J R Stat Soc Series B. 2007;69(2):243–268.  doi:10.1111/j.1467-9868.2007.00587.x.
  * [47] Bosse NI, Abbott S, Cori A, van Leeuwen E, Bracher J, Funk S. Transformation of forecasts for evaluating predictive performance in an epidemiological context.  medRxiv [Preprint]; 2023 [cited 2022 March 3].  Available from: https://doi.org/10.1101/2023.01.23.23284722. doi:10.1101/2023.01.23.23284722.
  * [48] Bosse N, Abbott S, Gruson H, Bracher J, Funk S. Scoringutils: Utilities for Scoring and Assessing Predictions; 2022.  R package version 1.0.1.  Available from: https://cran.r-project.org/web/packages/scoringutils.
  * [49] Funk S, Camacho A, Kucharski AJ, Lowe R, Eggo RM, Edmunds WJ.  Assessing the performance of real-time epidemic forecasts: A case study of Ebola in the Western Area region of Sierra Leone, 2014-15.  PLoS Comput Biol. 2019;15(2):e1006785.  doi:10.1371/journal.pcbi.1006785.
  * [50] Jordan A, Krueger F, Lerch S. ScoringRules: Scoring Rules for Parametric and Simulated Distribution Forecasts; 2022.  R package version 1.0.2.  Available from: https://cran.r-project.org/web/packages/scoringRules.
  * [51] McDonald DJ, Bien J, Green A, Hu AJ, DeFries N, Hyun S, et al.  Can auxiliary indicators improve COVID-19 forecasting and hotspot prediction?  Proc Natl Acad Sci U S A. 2021;118(51):e2111453118.  doi:10.1073/pnas.2111453118.
  * [52] R Core Team. R: A language and environment for statistical computing.  R Foundation for Statistical Computing; 2022.  Version 4.2.1.  Available from: https://www.R-project.org/.
  * [53] Brauner JM, Mindermann S, Sharma M, Johnston D, Salvatier J, Gavenčiak T, et al.  Inferring the effectiveness of government interventions against COVID-19.  Science. 2021;371(6531):eabd9338.  doi:10.1126/science.abd9338.
  * [54] Banholzer N, Weenen Ev, Lison A, Cenedese A, Seeliger A, Kratzwald B, et al.  Estimating the effects of non-pharmaceutical interventions on the number of new infections with COVID-19 during the first epidemic wave.  PLoS One. 2021;16(6):e0252827.  doi:10.1371/journal.pone.0252827.
  * [55] Cramer EY, Ray EL, Lopez VK, Bracher J, Brennen A, Castro Rivadeneira AJ, et al.  Evaluation of individual and ensemble probabilistic forecasts of COVID-19 mortality in the United States.  Proc Natl Acad Sci U S A. 2022;119(15):e2113561119.  doi:10.1073/pnas.2113561119.
  * [56] Sherratt K, Gruson H, Grah R, Johnson H, Niehus R, Prasse B, et al.. Predictive performance of multi-model ensemble forecasts of COVID-19 across European nations.  medRxiv [Preprint]; 2022 [cited 2022 March 3].  Available from: https://www.medrxiv.org/content/10.1101/2022.06.16.22276024v1. doi:10.1101/2022.06.16.22276024.
  * [57] Reinhart A, Brooks L, Jahja M, Rumack A, Tang J, Agrawal S, et al.  An open repository of real-time COVID-19 indicators.  Proc Natl Acad Sci U S A. 2021;118(51):e2111452118.  doi:10.1073/pnas.2111452118.
  * [58] Schwabe A, Persson J, Feuerriegel S.  Predicting COVID-19 spread from large-scale mobility data.  In: Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining; 2021. p. 3531–3539.  Available from: https://doi.org/10.1145/3447548.3467157. doi:10.1145/3447548.3467157.
  * [59] Rodríguez A, Tabassum A, Cui J, Xie J, Ho J, Agarwal P, et al.  DeepCOVID: An Operational Deep Learning-driven Framework for Explainable Real-time COVID-19 Forecasting.  In: Proceedings of the AAAI Conference on Artificial Intelligence. vol. 35; 2021. p. 15393–15400.  Available from: https://ojs.aaai.org/index.php/AAAI/article/view/17808. doi:10.1609/aaai.v35i17.17808.
  * [60] Reiner RC, Barber RM, Collins JK, Zheng P, Adolph C, Albright J, et al.  Modeling COVID-19 scenarios for the United States.  Nat Med. 2021;27(1):94–105.  doi:10.1038/s41591-020-1132-9.
  * [61] Noh J, Danuser G.  Estimation of the fraction of COVID-19 infected people in U.S. states and countries worldwide.  PLoS One. 2021;16(2):e0246772.  doi:10.1371/journal.pone.0246772.

## Acknowledgments

We thank Anne Cori for her support with the implementation of EpiEstim. We
thank Nikos Bosse for his suggestions on evaluating probabilistic forecasts.

## Ethics approval

Ethics approval was not required for this study.

## Competing interests

SF reports grants from the Swiss National Science Foundation outside of the
submitted work. All other authors declare no competing interests.

## Funding

NB and SF acknowledge funding from the Swiss National Science Foundation
(SNSF) as part of the Eccellenza grant 186932 on "Data-driven health
management". TM, SM, HJTU and SB acknowledge funding from the MRC Centre for
Global Infectious Disease Analysis (reference MR/R015600/1), jointly funded by
the UK Medical Research Council (MRC) and the UK Foreign, Commonwealth and
Development Office (FCDO), under the MRC/FCDO Concordat agreement and is also
part of the EDCTP2 programme supported by the European Union; and acknowledges
funding by Community Jameel. The funding bodies had no control over design,
conduct, data, analysis, review, reporting, or interpretation of the research
conducted. SB is funded by the National Institute for Health Research (NIHR)
Health Protection Research Unit in Modelling and Health Economics, a
partnership between the UK Health Security Agency, Imperial College London and
LSHTM (grant code NIHR200908). Disclaimer: “The views expressed are those of
the author(s) and not necessarily those of the NIHR, UK Health Security Agency
or the Department of Health and Social Care.” SB acknowledges financial
support from the Novo Nordisk Foundation via The Novo Nordisk Young
Investigator Award (NNF20OC0059309) which also supports SM. SB acknowledges
financial support from the Danish National Research Foundation via a chair
grant. SB acknowledges financial support from The Eric and Wendy Schmidt Fund
For Strategic Innovation via the Schmidt Polymath Award (G-22-63345).

## Author Contributions

NB, TM, SM, and SB conceptualized the study. NB collected and analyzed the
data. NB, TM, SM, and SB contributed to methodology. All authors contributed
to interpretation and visualization of data and results. NB wrote the original
draft. All authors contributed to writing – review & editing.

## Supporting information

S1 Appendix. Supplementary material. Fig S1. Labeling of the epidemic
trajectories. Fig S2. Forecasting performance for different modeling choices.
Fig S3. Probabilistic forecasts for Arizona by model. Fig S4. Probabilistic
forecasts for Illinois by model. Fig S5. Probabilistic forecasts for Maryland
by model. Fig S6. Probabilistic forecasts for New Jersey by model. Fig S7.
Probabilistic forecasts for New York by model. Fig S8. Evaluation of
forecasting performance by US state. Fig S9. Evaluation of forecasting
performance by epidemic phase.

S1 Appendix

Fig S1: Labeling of the epidemic trajectories. Epidemic trajectories in US
states manually labeled according to five phases: exponential increase,
subexponential increase, plateau, subexponential decline, exponential decline.
Fig S2: Forecasting performance for different modeling choices. CRPS (left
average scores, right boxplot of scores) for the 1- and 2-week ahead forecast
by model specification. (a) EpiEstim. (b) EpiNow2. (c) SARIMA. (d) Prophet.
(e) GP. Fig S3: Probabilistic forecasts for Arizona by model. Probabilistic
forecast (posterior mean as red line, 95%-PI as shaded area) for Arizona by
each model. (a) 1-week ahead forecast. (b) 2-week ahead forecast. Observed
incidence is shown with black lines. Fig S4: Probabilistic forecasts for
Illinois by model. Probabilistic forecast (posterior mean as red line, 95%-PI
as shaded area) for Illinois by each model. (a) 1-week ahead forecast. (b)
2-week ahead forecast. Observed incidence is shown with black lines. Fig S5:
Probabilistic forecasts for Maryland by model. Probabilistic forecast
(posterior mean as red line, 95%-PI as shaded area) for Maryland by each
model. (a) 1-week ahead forecast. (b) 2-week ahead forecast. Observed
incidence is shown with black lines. Fig S6: Probabilistic forecasts for New
Jersey by model. Probabilistic forecast (posterior mean as red line, 95%-PI as
shaded area) for New Jersey by each model. (a) 1-week ahead forecast. (b)
2-week ahead forecast. Observed incidence is shown with black lines. Fig S7:
Probabilistic forecasts for New York by model. Probabilistic forecast
(posterior mean as red line, 95%-PI as shaded area) for New York by each
model. (a) 1-week ahead forecast. (b) 2-week ahead forecast. Observed
incidence is shown with black lines. Fig S8: Evaluation of forecasting
performance by US state. Median CRPS by US state for the 1- and 2-week ahead
forecast by model. Fig S9: Evaluation of forecasting performance by epidemic
phase. Median CRPS by different epidemic phases (see labeling in Figure S1)
for the 1- and 2-week ahead forecast by model.