[a]Sofie Martins

# Finite-Size Effects of the HVP Contribution to the Muon $g-2$ with C⋆
Boundary Conditions

Agostino Patella

###### Abstract

The muon $g-2$ is a compelling quantity due to the current standing tensions
among the experimental average, data-driven theoretical results, and lattice
results. Matching the final target accuracy of the experiments at Fermilab and
J-PARC will constitute a major challenge for the lattice community in the
coming years. For this reason, it is worthwhile to consider different options
to keep the systematic errors under control. In this proceedings, we discuss
finite-volume effects of the leading Hadron Vacuum Polarization (HVP)
contribution to the muon $g-2$ in the presence of C⋆ boundary conditions. When
considering isospin-breaking corrections to the HVP, C⋆ boundary conditions
provide a possible consistent formulation of $\mathrm{QCD+QED}$ in finite
volume. Even though these boundary conditions can be avoided in the
calculation of the leading HVP contribution, we find the interesting result
that they remove the leading exponential finite-volume correction. In
practice, compared to the periodic case, C⋆ boundary conditions cut the
finite-size effects in half on a lattice of physical size $M_{\pi}L=4$ and by
a factor of almost ten for $M_{\pi}L=8$. We discuss the origin of this
reduction and implications for computational efficiency.

##  1 Motivation

In this proceedings, we discuss finite-volume effects of the HVP contribution
to the muon $g-2$ in isosymmetric QCD with C⋆ boundary conditions [1, 2, 3,
4]. A priori, there is no reason to choose these boundary conditions for QCD
observables. Nevertheless, it is interesting to consider this setup since we
are interested in using it to calculate eventually also isospin-breaking
corrections to the HVP. As discussed in [5], C⋆ boundary conditions provide
the only finite-volume formulation of QED which preserves locality, gauge
invariance and translational invariance [5]. First numerical explorations in
this direction have been presented in [6, 7, 8]. We recall that locality can
also be preserved by adding a mass term to the photon while breaking gauge
invariance in a controlled way [9, 10].

Interestingly, we find that C⋆ boundary conditions reduce the finite-volume
correction to the HVP contribution to the muon $g-2$. In fact, these are order
$e^{-M_{\pi}L}$ and $e^{-\sqrt{2}M_{\pi}L}$ with periodic and C⋆ boundary
conditions, respectively. While the proof of this statement will be given in a
future publication, we sketch here the main ideas and the most notable
differences with the periodic case, analyzed in [11, 12]. We also provide an
estimate of the finite-size corrections with C⋆ boundary conditions using a
phenomenological model for the pion form factor.

##  2 C⋆ Boundary Conditions

### Definition.

[5] introduces the simulation of charged hadrons on the lattice using C⋆
boundary conditions. Quarks behave as in the following

$\Psi_{f}(x+L\hat{e}_{i})=\Psi_{f}(x)=C^{-1}\bar{\Psi}_{f}^{\mathrm{T}}(x)\,,$
(1)
$\bar{\Psi}_{f}(x+L\hat{e}_{i})=\bar{\Psi}_{f}(x)=-\Psi_{f}^{\mathrm{T}}(x)C\,,$
(2)

where $L$ is the lattice extent in space. Fields that are eigenstates of the
charge conjugation operator behave as periodic or antiperiodic fields
depending on their C-parity. For example, the electromagnetic current has
negative C-parity and is an antiperiodic field in C⋆ boundary conditions:

$j_{\mu}(x+L\hat{e}_{i})=j_{\mu}^{\mathrm{c}}(x)=-j_{\mu}(x)\,.$ (3)

As a consequence, the allowed momenta for the electromagnetic current are
given by

$\Pi_{-}=\left\\{\dfrac{\pi}{L}(2\bm{n}+\bar{\bm{n}})|\bm{n}\in\mathbb{Z}^{3},\bar{\bm{n}}=(1,1,1)\right\\}\,$
(4)

see [5]. This property has consequences for the definition of a finite-volume
estimator for the leading-order HVP. While in the periodic case, one simply
considers the zero-momentum two-point function of the electromagnetic current,
this is not an option in the case of C⋆ boundary conditions, since the
electromagnetic current does not have a zero-momentum component.

### Finite-volume estimator.

Irrespectively of the boundary conditions, a finite-volume estimator of the
leading-order HVP contribution to the muon $g-2$ can be defined as

$a_{\mu}^{\mathrm{LO,HVP}}=\int_{0}^{T/2}\mathrm{d}x_{0}\,\mathcal{K}(x_{0})G(x_{0}|T,L)\,,$
(5)

where $\mathcal{K}(x_{0})$ is the kernel defined in [13], and we define

$G(x_{0}|T,L)=-\dfrac{1}{3}\int_{V_{L}}\langle
j_{k}(x)j_{k}(0)\rangle_{T,L}\mathrm{d}^{3}x\,.$ (6)

The integral of the two-point function is taken over a cube $V_{L}$ with an
edge length equal to $L$. In the case of periodic boundary conditions, since
the electromagnetic current is periodic in space, the integration domain
$V_{L}$ can be freely translated without changing the value of the integral.
In particular, the integral corresponds to a zero-momentum projection. In the
case of C⋆ boundary conditions, since the electromagnetic current is
antiperiodic, translating the integration domain $V_{L}$ will change the value
of the integral, yielding a different finite-volume estimator. One can show
that the symmetric choice

$V_{L}=\left(-\dfrac{L}{2},\dfrac{L}{2}\right)^{3}$ (7)

has the correct infinite-volume limit and yields finite-volume effects that
are smaller with respect to the periodic case. Notice that, in the case of C⋆
boundary conditions, the integral over $V_{L}$ does not correspond to the
projection over a single momentum: states propagating in between the two
currents are superimpositions of all allowed momenta.

In principle, different finite-volume estimators can be constructed, e.g., by
projecting the two-point function on one of the momenta allowed by the
boundary conditions. For instance one could choose
$\bm{p}=\frac{\pi}{L}(1,1,1).$ However, we found that this generates finite-
volume corrections that vanish as inverse powers of $L$. For this reason, this
choice will not be considered further.

##  3 Finite-Size Effects in the Isospin-Symmetric Limit

To derive an analytical formula for the finite-size effects for C⋆ boundary
conditions, we can use the method developed in [12] and apply it with modified
boundary conditions. This calculation is non-trivial and beyond the scope of
this proceedings. However, it is possible to examine the structure of the
effects arising for C⋆ boundary conditions and compare them to periodic
boundary conditions.

### Necessary Adjustments.

[12] derives the finite-volume effects by expanding to all orders in an
effective theory of pions. While in periodic boundary conditions, all pions
are periodic fields, in C⋆ boundary conditions, their behavior is modified as
in the following

$\pi^{0}(x+L\hat{e}_{i})=\pi^{0}(x)\,,$ (8)
$\pi^{q}(x+L\hat{e}_{i})=\pi^{-q}(x)\,.$ (9)

As a result, in contrast to periodic boundary conditions, we need to
distinguish between contributions associated with different pion charges in
the final formula .

### Analytic Results.

Employing these adjustments, we find the following analytic expression

$\begin{split}\Delta G_{L}(x_{0})=&-\sum_{\bm{n}\neq\bm{0}}\sum_{q=\\{0,\pm
1\\}}\chi_{q,\bm{n}}\int\dfrac{\mathrm{d}p_{3}}{2\pi}\dfrac{e^{-|\bm{n}|L\sqrt{M_{\pi}^{2}+p_{3}^{2}}}}{24\pi|\bm{n}|L}\int\dfrac{\mathrm{d}k_{3}}{2\pi}\>\mathrm{cos}(k_{3}x_{0})\mathrm{Re}T^{q}(-k_{3}^{2},-p_{3}k_{3})\\\
&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\mathcal{O}(e^{-\sqrt{2+\sqrt{3}}M_{\pi}L})\end{split}$
(10)

where

$\chi_{q,\bm{n}}=\dfrac{1+(-1)^{q\sum_{i}n_{i}}}{2}\,,$ (11)

and $T^{q}(-k_{3}^{2},-p_{3}k_{3})$ is the forward Compton scattering
amplitude of a pion with charge $q$ and momentum $(0,0,p_{3})$ against an off-
shell photon with Minkowskian four-momentum $(0,0,0,k_{3})$, traced over the
Lorenz-index of the photon. The adjustments due to the modified boundary
conditions translate into this formula in the following way: Firstly, there is
an explicit sum over different charges and, therefore, contributions that are
specifically associated with the differently charged pions. For this, we need
to calculate the Compton scattering amplitude of differently charged pions,
here denoted $T^{q}$ instead of an implicitly summed pion Compton scattering
amplitude $T$.

Secondly, the factor $\chi$ needs to be specifically introduced for C⋆
boundary conditions. For the contribution through the uncharged pion $\chi=1$
for all $\bm{n}$, identically to periodic boundary conditions.

For the charged contributions, this factor can be either zero or one,
depending on the vector $\bm{n}$. This vector can be interpreted as counting
the number of times an internal pion line wraps around the torus, which is
defined as $\sum_{i=1}^{3}n_{i}$. For $q\neq 0$, any odd wrapping number will
yield $\chi=0$. As a result, contributions with odd $\bm{n}$ are filtered out,
and only contributions for $q=0$ at the odd orders remain. This affects most
notably the leading order.

### Spectral Decomposition.

These uncharged contributions now turn out to be small already. This is
because a spectral decomposition as in [12] yields three contributions, a
vacuum contribution, a one-particle pole contribution and remaining regular
contributions.

$T^{q}=T^{q}_{\mathrm{vac}}+T^{q}_{\mathrm{pole}}+T^{q}_{\mathrm{reg}}$ (12)

As in [12], we can assess all regular contributions, which is the regular part
of the one-particle contribution and the multi-particle and higher-mass terms
through chiral perturbation theory and find them to be small compared to all
other effects. They are included in our numerical estimates.

The most significant part is the pole contribution through one-particle
intermediate states. This term is proportional to the infinite-volume pion
form factor and contributes to the finite-volume effects for the charged
cases, which are removed through $\chi$ at odd orders. The remaining uncharged
contributions are small because the form factor is zero for $q=0$.

$T_{1\pi,\mathrm{pole}}^{0}\propto|\bra{\pi^{0}}J_{\mu}(0)\ket{\pi^{0}}|^{2}=0$
(13)

Further, the vacuum contribution is proportional to the matrix element

$T_{\mathrm{vac}}\propto\bra{\pi^{q}}J_{\mu}\ket{\Omega}$ (14)

where $\ket{\Omega}$ denotes the vacuum. This is zero for charged pions thanks
to charge conservation, and for the neutral pion thanks to C-parity
conservation.

### General Structure.

As for the periodic case discussed in [11] the effects exhibit exponential
suppression with pion mass $M_{\pi}$ and lattice extent $L$, now with odd
orders removed. Different orders of effects are generated by the sum over
$\bm{n}$, generating terms of the structure

$\mathcal{O}(e^{-\sqrt{2}M_{\pi}L})+\mathcal{O}(e^{-2M_{\pi}L})+\mathcal{O}(e^{-\sqrt{6}M_{\pi}L})+\ldots$
(15)

followed by mixed finite-time and finite-volume effects

$\mathcal{O}(e^{-M_{\pi}T})+\ldots+\mathcal{O}(e^{-M_{\pi}\sqrt{L^{2}+T^{2}}})+\ldots$
(16)

which for a setup of $T=2L$ vanish like $e^{-2M_{\pi}L}$. Finally, there are
contributions through higher-mass hadrons

$\mathcal{O}(e^{-\sqrt{2}M_{K}L})+\mathcal{O}(e^{-2M_{K}L})+\ldots$ (17)

which are also negligible because the next-higher-mass hadron is the kaon
which has $M_{K}>2M_{\pi}$. Terms in 16 and 17 are both neglected together
with other contributions in the final formula for finite-volume effects at the
order $\mathcal{O}(e^{-\sqrt{2+\sqrt{3}}M_{\pi}L})$.

### Periodic Boundary Conditions.

The main difference to the periodic boundary conditions from [12] is the
necessity to introduce the factor $\chi$. One can recover the analytic result
for periodic boundary conditions by setting $\chi=1$ in eq. 10.

##  4 Results and Conclusion

Figure 1: Finite-Volume Effects are cut in half for a physical lattice size of
$M_{\pi}L=4$ when using C⋆ boundary conditions compared to periodic boundary
conditions. For a large lattice of size $M_{\pi}L=8$ the effects are cut by a
factor of almost ten. We are using a pion mass at the physical point of
$M_{\pi}=137\,\mathrm{MeV}$. The finite-volume effects contribute at permil
for C⋆ boundary conditions for a lattice of physical size $M_{\pi}L\approx 6$
while for periodic boundary conditions $M_{\pi}L\approx 8$ is necessary

(a) Absolute Finite-Volume Effects to the HVP $M_{\pi}L$ | C⋆ BC | PBC[12]
---|---|---
4 | 9.74(1.6) | 22.4(3.1)
5 | 3.25(0.23) | 10.0(0.4)
6 | 1.027(0.034) | 4.42(0.06)
7 | 0.311(0.005) | 1.924(0.009)
8 | 0.0909(0.0008) | 0.826(0.001)

(b) Relative Finite-Volume Effects to the HVP $M_{\pi}L$ | C⋆ BC | PBC[12]
---|---|---
4 | 1.39 | 3.20
5 | 0.464 | 1.43
6 | 0.147 | 0.631
7 | 0.0444 | 0.275
8 | 0.0130 | 0.118

1%0.1%0.01%$0$$10$$20$$30$$40$Relative Contribution of FV EffectsEstimated
Cost FactorC⋆ BCPBC Figure 2: Comparison of Finite-Size Effects Reached with
Different Boundary Conditions

### Comparison of Finite-Volume Effects in Different Lattice Sizes.

We can evaluate these analytic formulae numerically analogously to the
approach taken in [12]. Table 1 shows the absolute contributions of finite-
volume effects to the HVP, comparing periodic boundary conditions (PBC) and C⋆
boundary conditions (C⋆ BC) for different physical lattice sizes. For a small
lattice of $M_{\pi}L=4$, the finite-volume effects are cut by more than half
compared to periodic boundary conditions. For a larger lattice of
$M_{\pi}L=8$, this effect is even stronger: The effects are reduced by almost
a factor of ten. From comparing the relative contributions in table 1, one can
see permil finite-size effects are reached for periodic boundary conditions
for $M_{\pi}L=8$ while for C⋆ boundary conditions, a lattice of physical size
$M_{\pi}L=6$ is already sufficient.

### Estimating Performance Differences.

When comparing the periodic and C⋆ setup, one needs to take into account the
fact that simulations with C⋆ boundary conditions are generally more
expensive. This is partly due to the fact that the lattice needs to be
effectively doubled, and partly to the fact that the HMC algorithm for the
light quarks needs to be replaced by the RHMC. We estimate that isosymmetric
QCD simulations with C⋆ boundary conditions are about 2.3 times more expensive
than simulations with periodic boundary conditions. Considering the
computational cost, we can compare the merit of the different boundary
conditions at different target finite-volume effects. Figure 2 shows an
estimate of the computational time necessary to achieve a certain size in the
finite-size effects comparing the two choices for boundary conditions. At a
fixed relative size of the finite-volume effects, starting from 1%, we expect
C⋆ boundary conditions to outperform periodic boundary conditions
consistently; however, for smaller target finite-volume effects, this
improvement is marginal. Since we are targeting permil precision for the HVP,
C⋆ boundary conditions are plausibly the better choice for evaluating the
leading-order contribution.

### Outlook.

In this proceedings we have described finite-volume corrections to the HVP
contribution to the muon $g-2$, calculated in isosymmetric QCD with C⋆
boundary conditions. The proof of the formulae presented in this proceedings
will be given in a future publication. The assessment of the finite-size
effects in the isospin-breaking corrections to the HVP will be the subject of
future work.

### Acknowledgements.

S.M. received funding from the European Union’s Horizon 2020 research and
innovation program under the Marie Skłodowska-Curie grant agreement № 813942.

## References

  * Kronfeld and Wiese [1991] A. S. Kronfeld and U. J. Wiese, Nucl. Phys. B 357, 521 (1991).
  * Kronfeld and Wiese [1993] A. S. Kronfeld and U. J. Wiese, Nucl. Phys. B 401, 190 (1993), arXiv:hep-lat/9210008 .
  * Polley and Wiese [1991] L. Polley and U. J. Wiese, Nucl. Phys. B 356, 629 (1991).
  * Wiese [1992] U. J. Wiese, Nucl. Phys. B 375, 45 (1992).
  * Lucini  _et al._ [2016] B. Lucini, A. Patella, A. Ramos, and N. Tantalo, JHEP 02, 076 (2016), arXiv:1509.01636 [hep-th] .
  * Tavella [2022] P. Tavella, “Strange and charm contribution to the hvp from c* boundary conditions,”  (2022), talk at The 39th International Symposium on Lattice Field Theory.
  * Gruber [2022] R. Gruber, “A first look at the hvp from qcd and qcd+qed with c* boundary conditions,”  (2022), talk at The 39th International Symposium on Lattice Field Theory.
  * Altherr [2022] A. Altherr, “A first look at the hvp from qcd and qcd+qed with c* boundary conditions ii,”  (2022), talk at The 39th International Symposium on Lattice Field Theory.
  * Clark  _et al._ [2022] M. A. Clark, M. Della Morte, Z. Hall, B. Hörz, A. Nicholson, A. Shindler, J. T. Tsang, A. Walker-Loud, and H. Yan, PoS LATTICE2021, 281 (2022), arXiv:2201.03251 [hep-lat] .
  * Endres  _et al._ [2016] M. G. Endres, A. Shindler, B. C. Tiburzi, and A. Walker-Loud, Phys. Rev. Lett. 117, 072002 (2016), arXiv:1507.08916 [hep-lat] .
  * Hansen and Patella [2019] M. T. Hansen and A. Patella, Phys. Rev. Lett. 123, 172001 (2019), arXiv:1904.10010 [hep-lat] .
  * Hansen and Patella [2020] M. T. Hansen and A. Patella, JHEP 10, 029 (2020), arXiv:2004.03935 [hep-lat] .
  * Bernecker and Meyer [2011] D. Bernecker and H. B. Meyer, Eur. Phys. J. A 47, 148 (2011), arXiv:1107.4388 [hep-lat] .


# Spatial forecast postprocessing: The Max-and-Smooth approach

Stefan Siegert111Corresponding author, Email<EMAIL_ADDRESS>
Ben Hooper, Joshua Lovegrove, Tyler Thomson,
Birgir Hrafnkelsson

(Draft version: )

## Abstract

Numerical weather forecasts can exhibit systematic errors due to simplifying
model assumptions and computational approximations. Statistical postprocessing
is a statistical approach to correcting such biases. A statistical
postprocessing model takes input data from a numerical forecast model, and
outputs a parametric predictive distribution of a real-world observation, with
model parameters learned from past forecast-observation pairs. In this paper
we develop and discuss methods for postprocessing of gridded data. We show
that estimates of postprocessing parameters on a spatial grid can be improved
by Bayesian hierarchical modelling with spatial priors. We use the “Max-and-
Smooth” approach (Hrafnkelsson et al., 2021) to approximate a fully Bayesian
inference in two steps. First we calculate maximum-likelihood estimates (MLEs)
of postprocessing parameters at individual grid points. Second we smooth the
MLEs using a measurement error model with a spatial prior. Our approach
provides the theoretical basis for the parameter smoothing approach by Kharin
et al. (2017), and simplifies and generalises the Bayesian hierarchical
modelling approach by Möller et al. (2015). A new derivation of Max-and-Smooth
is presented. The method is applicable to arbitrary postprocessing models, as
illustrated on Model Output Statistics, Logistic Regression, and
Nonhomogeneous Gaussian Regression. We report consistent improvements in
forecast accuracy, calibration, and probabilistic skill in postprocessing of
temperature and precipitation forecasts.

##  1 Introduction

Forecasts generated by numerical simulators are prone to systematic biases
that can be corrected in a postprocessing step (Glahn and Lowry, 1972;
Vannitsem et al., 2018). A common example is to subtract a constant bias from
a temperature forecast that is known to be systematically too warm.
Statistical postprocessing techniques use a data set of past forecasts, and
their known verifying observations, to fit a parametric statistical model that
takes a forecast as input and outputs a point estimate or a probability
distribution for the corresponding future observation. Due to the dynamic
nature of the forecast model and underlying physical system, the statistical
behavior of forecast errors can change over space and time. Postprocessing
models should be able to adapt to such local differences, e.g., by using
different statistical models in different locations. The research presented in
this paper addresses the problem of efficiently using spatial forecast and
observation data to fit postprocessing models with spatially adaptive
parameters.

We differentiate between spatial postprocessing and local postprocessing.
Local postprocessing approaches fit statistical postprocessing models
independently at each location, using only data from that location, without
accounting for any spatial correlations of data or model parameters. In
spatial postprocessing statistical models are used that can account for
spatial correlations; see Schefzik and Möller (2018) for a review. Spatial
postprocessing can be applied either at the data level, or the parameter
level, or both. Spatial postprocessing at the data level, i.e., applied at the
level of the predicted meteorological variable, aims to produce spatially
coherent meteorological fields. Data-level methods include the Schaake shuffle
(Clark et al., 2004), ensemble copula coupling (Schefzik et al., 2013), and
Gaussian process models (Berrocal et al., 2008). Spatial postprocessing at the
parameter level aims to improve estimates of local postprocessing parameters
by exploiting similarity of parameters at nearby locations. Approaches include
postprocessing of spatially aggregated data (van Straaten et al., 2020),
spatial smoothing of postprocessing parameters (Kharin et al., 2017; Lovegrove
and Siegert, 2022), or Bayesian hierarchical modelling with a spatially
correlated prior on postprocessing parameters (Möller et al., 2015).
Postprocessing at parameter and data level can be combined to improve both,
estimates of local model parameter and the spatial coherence of output fields
(Möller et al., 2015).

This paper is concerned with spatial postprocessing at the parameter level. We
present a simple method improve local estimates of postprocessing parameters
by borrowing strength from data at nearby locations. Formally the method
yields the same results as the Max-and-Smooth approach by Hrafnkelsson et al.
(2021), but the derivation of the method differs slightly from the original
proof. The argument starts from a Bayesian hierarchical model for the
meteorological observations, with a spatially correlated prior distribution
used for the postprocessing parameters. In that way our method is similar to
the approach by Möller et al. (2015). We show that the posterior distribution
of postprocessing parameters can be approximated by first calculating their
local maximum likelihood estimates (MLEs), and secondly, spatially smoothing
the MLEs. In that way our method is similar to the parameter smoothing method
proposed by Kharin et al. (2017). The “Max-and-Smooth” approach at the same
time unifies the approaches of Möller et al. (2015) and Kharin et al. (2017),
it simplifies and speeds up the Bayesian inference used by Möller et al.
(2015), and provides a more solid mathematical foundation for the ad-hoc
smoothing approach of Kharin et al. (2017). In our method, the degree of
spatial smoothing is partly controlled by the estimation uncertainty of the
local MLEs, which allows us to learn the smoothing hyperparameters from the
data, which is another advantage compared to the approach of Kharin et al.
(2017).

Temperature and precipitation data used to illustrate our methods are detailed
in Sec. 2.1. We illustrate the Max-and-Smooth approach on postprocessing
parameters from Model Output Statistics (MOS), Logistic Regression (LR), and
Nonhomogeneous Gaussian Regression (NGR), which are briefly summarised in Sec.
2.2. The main result of the paper is the Bayesian Hierarchical modelling
approach to spatial postprocessing and its formulation as a spatial smoothing
applied to Maximum likelihood estimates of postprocessing parameters, which is
detailed in Sec. 2.3. Further methodological details are summarised in Sec.
2.4 on choice of the spatial prior, in Sec. 2.5 on inference of
hyperparameters, Sec. 2.6 on the R-INLA software package, and Sec. 2.7 on
evaluation methods and metrics used to assess the performance of the proposed
method. Sec. 3 summarises performance improvements due to spatial
postprocessing of temperature and precipitation forecasts with LR, MOS and
NGR, and Sec. 4 concludes with a summary and discussion.

##  2 Data and Methods

###  2.1 Data

To illustrate and compare the methods we use hindcast data of different
meteorological variables and at different forecast lead times. Forecast data
were downloaded from the ECMWF S2S data base (Vitart et al., 2017). Forecasts
were generated by ECMWF’s seasonal forecast model which has 10 ensemble
members perturbed from 1 control forecast. Reanalysis data from the ERA5
reanalysis project (Hersbach et al., 2020) were used as pseudo-observations.
All forecast and reanalysis data are available on a $1.5^{\circ}$ latitude-
longitude grid on the rectangular domain between longitudes $15^{\circ}W$ and
$30^{\circ}E$ and latitudes $36^{\circ}N$ and $69^{\circ}N$. The API requests
to download reanalysis and reforecast data are given in the Supplementary
Material.

We use 2-metre temperature ensemble forecasts to illustrate and compare
performance metrics when postprocessing real-valued quantities. We use
exceedance of precipitation thresholds to illustrate the methods applied to
postprocessing probabilities of binary events. All forecasts were initialised
at 00:00UTC on 14 April each year from 2002 to 2021, and forecast lead times 1
day up to 10 days.

###  2.2 Postprocessing methods

We assume hindcast data at a given lead time are available as ensemble
forecasts with $K$ members on a regular latitude-longitude grid with grid
point indexed by $s=1,\dots,S$, and for verification date indexed by
$t=1,\dots,T$. Since all analyses are carried out independently for each
forecast lead time, the lead time index is suppressed. We denote the $k$-th
ensemble member at grid point $s$ and verification date $t$ by $f_{k,s,t}$,
and summarise the ensemble by the ensemble mean

$m_{s,t}=\frac{1}{K}\sum_{k=1}^{K}f_{k,s,t}$

and the ensemble variance

$v_{s,t}=\frac{1}{K-1}\sum_{k=1}^{K}(f_{k,s,t}-m_{s,t})^{2}.$

Model output statistics (MOS, Glahn and Lowry, 1972) was one of the earliest
postprocessing methods proposed to correct systematic biases of numerical
models, and to predict variables that are not explicitly simulated by the
numerical model. MOS is essentially a least-squares linear regression, where
the verifying observation $y_{s,t}$ is modelled as a linear transformation of
the ensemble mean forecast $m_{s,t}$, plus an independent normally distributed
error with zero mean and time-constant variance:

$y_{s,t}\sim N(\mu^{MOS}_{s,t},(\sigma^{MOS}_{s})^{2})$ (1)

where

$\displaystyle\mu^{MOS}_{s,t}$
$\displaystyle=\alpha^{MOS}_{s}+\beta^{MOS}_{s}(m_{s,t}-\bar{m}_{s}),$ (2)
$\displaystyle\log\left(\sigma^{MOS}_{s}\right)^{2}$
$\displaystyle=\tau^{MOS}_{s},$ (3)

and where $\bar{m}_{s}=\frac{1}{T}\sum_{t=1}^{T}m_{s,t}$ is the time-averaged
ensemble mean. We parameterise the model in terms of the log-variance
$\tau^{MOS}_{s}$ to have a real-valued parameter for which we can use a normal
distribution as a spatial prior for Bayesian inference. Centering the
covariates by subtracting $\bar{m}_{s}$ diagonalises the information matrix
(below) and hence makes MOS parameter estimates asymptotically independent.

Let $\boldsymbol{y}_{s}$ denote the vector of observed data at location $s$.
The log-likelihood function of the MOS model is given by

$\displaystyle\log p(\boldsymbol{y}_{s}|\boldsymbol{\theta}^{MOS}_{s})=\log
p(\boldsymbol{y}_{s}|\alpha^{MOS}_{s},\beta^{MOS}_{s},\tau^{MOS}_{s})$ (4)
$\displaystyle=-\frac{T}{2}\log
2\pi-\frac{T}{2}\tau^{MOS}_{s}-\frac{1}{2}e^{-\tau^{MOS}_{s}}\sum_{t=1}^{T}\left\\{y_{s,t}-(\alpha^{MOS}_{s}+\beta^{MOS}_{s}(m_{s,t}-\bar{m}_{s}))\right\\}^{2}.$

The maximum likelihood estimators (MLEs) at grid point $s$ are then given by

$\displaystyle\hat{\alpha}^{MOS}_{s}$
$\displaystyle=\bar{y}_{s}=\frac{1}{T}\sum_{t=1}^{T}y_{s,t},$ (5)
$\displaystyle\hat{\beta}^{MOS}_{s}$
$\displaystyle=\frac{\sum_{t=1}^{T}(y_{s,t}-\bar{y}_{s})(m_{s,t}-\bar{m}_{s})}{\sum_{t=1}^{T}(m_{s,t}-\bar{m}_{s})^{2}}\text{,
and}$ $\displaystyle\hat{\tau}^{MOS}_{s}$
$\displaystyle=\log\left\\{\frac{1}{T}\sum_{t=1}^{T}\left[y_{s,t}-(\hat{\alpha}^{MOS}_{s}+\hat{\beta}^{MOS}_{s}(m_{s,t}-\bar{m}_{s}))\right]^{2}\right\\}.$

We estimate the precision of the sampling distribution of the MLEs by the
observed information matrix (Efron and Hinkley, 1978) which is given by

$\hat{J}_{s}=-\frac{\partial^{2}\log
p(\boldsymbol{y}_{s}|\boldsymbol{\theta}_{s})}{\partial\boldsymbol{\theta}_{s}\partial\boldsymbol{\theta}_{s}^{\prime}}\Bigg{|}_{\boldsymbol{\theta}_{s}=\hat{\boldsymbol{\theta}}_{s}}=\begin{pmatrix}Te^{-\hat{\tau}_{s}}&0&0\\\
0&\sum_{t=1}^{T}(m_{s,t}-\bar{m}_{t})^{2}e^{-\hat{\tau}_{s}}&0\\\
0&0&\frac{T}{2}\end{pmatrix}.$ (6)

Logistic regression (Wilks, 2011, ch. 7) can be used to calculate probability
of precipitation from deterministic model output of total precipitation
amounts. We use simple logistic regression to predict the event that the daily
precipitation amount exceeds $2.5$mm, using the ensemble mean total
precipitation as a single covariate. Let the binary prediction target
$y_{s,t}$ be equal to $1$ if precipitation amount at location $s$ and
verification time $t$ exceeds $2.5$mm, and $y_{s,t}=0$ otherwise. Define the
forecast probability $p^{LR}_{s,t}=P(y_{s,t}=1)$, which is calculated from the
ensemble mean precipitation amount forecast $m_{s,t}$ by

$p^{LR}_{s,t}=\left(1+e^{-(\alpha^{LR}_{s}+\beta^{LR}_{s}m_{s,t})}\right)^{-1}.$
(7)

With $\boldsymbol{\theta}^{LR}_{s}=(\alpha^{LR}_{s},\beta^{LR}_{s})$, the log-
likelihood is given by

$\displaystyle\log p(\boldsymbol{y}_{s}|\boldsymbol{\theta}^{LR}_{s})$ (8)
$\displaystyle=-\sum_{t=1}^{T}\left[(1-y_{s,t})(\alpha^{LR}_{s}+\beta^{LR}_{s}m_{s,t})+\log\left(1+e^{-(\alpha^{LR}_{s}+\beta^{LR}_{s}m_{s,t})}\right)\right]$

There are no closed form solutions for the MLEs
$\hat{\boldsymbol{\theta}}^{LR}_{s}=(\hat{\alpha}_{s}^{LR},\hat{\beta}^{LR}_{s})$,
so they are obtained by numerical optimsation using the R function optim. The
observed information of the logistic regression model at location $s$ is given
by

$\hat{J}_{s}=\sum_{t=1}^{T}\hat{p}_{s,t}^{LR}(1-\hat{p}_{s,t}^{LR})\begin{pmatrix}1&m_{s,t}\\\
m_{s,t}&m_{s,t}^{2}\end{pmatrix}$ (9)

where $\hat{p}_{s,t}^{LR}$ is Eq. 7 evaluated at
$\hat{\boldsymbol{\theta}}_{s}^{LR}$.

Non-homogeneous Gaussian regression (NGR, Gneiting et al., 2005) is similar to
MOS but allows for non-constant forecast variances. NGR accounts for the fact
that in addition to the ensemble mean, the ensemble variance can be
informative for the forecast distribution. In NGR the real-valued observations
are normally distributed

$y_{s,t}\sim N(\mu^{NGR}_{s,t},(\sigma^{NGR}_{s,t})^{2})$ (10)

where the forecast mean $\mu^{NGR}_{s,t}$ depends linearly on the the ensemble
mean $m_{s,t}$ (as in MOS), and the forecast variance
$(\sigma^{NGR}_{s,t})^{2}$ depends linearly on the ensemble variance $v_{s,t}$
through

$\displaystyle\mu^{NGR}_{s,t}$ $\displaystyle=\alpha_{s}+\beta_{s}m_{s,t}$
(11) $\displaystyle\sigma^{NGR}_{s,t}$
$\displaystyle=(e^{\gamma_{s}}+e^{\delta_{s}}v_{s,t})^{1/2}.$

The log-likelihood follows from Eq. 10 and Eq. 11. The MLEs
$\hat{\boldsymbol{\theta}}^{NGR}_{s}=(\hat{\alpha}^{NGR}_{s},\hat{\beta}^{NGR}_{s},\hat{\gamma}^{NGR}_{s},\hat{\delta}^{NGR}_{s})^{\prime}$
are not available in closed form, and so we obtain them via numerical
optimisation. The observed information can be calculated analytically, but the
expressions are complicated, and we found that numerical approximations
(returned by optim) are accurate. Note that our implementation of NGR differs
from that of Gneiting et al. (2005), who used minimum CRPS estimation rather
than MLE, and parameterised the variance as $\gamma+\delta^{2}v_{s,t}$.

We are working with samples of $T=20$ when fitting local MLEs of
postprocessing parameters. Such small sample sizes can lead to poor
convergence and “exploding” parameter estimates in numerical optimisation. To
help the numerical optimiser and somewhat regularise the parameter estimates
we added $-0.5\times 10^{-4}[(\alpha^{LR})^{2}+(\beta^{LR})^{2}]$ to the LR
log-likelihood, and, similarly $-0.5\times
10^{-4}[(\alpha^{NGR})^{2}+(\beta^{NGR})^{2}+(\gamma^{NGR})^{2}+(\delta^{NGR})^{2}]$
to the NGR log-likelihood during numerical optimisation. This additive term
can be interpreted as an $L_{2}$ (ridge) penalty that shrinks parameter
estimates towards zero (Hastie et al., 2009), or equivalently as independent
zero-mean normal prior distributions with standard deviation 100 applied to
all parameters. While we will refer to parameter estimates as MLEs, strictly
speaking they should be called PMLEs (penalised MLEs) or MAP (maximum
a-posteriori) estimates. The effect of the penalty term is negligible in cases
where we optimiser converges without the penalty, and it leads to successful
convergence to “reasonable” parameter values when the optimiser fails to
converge without the penalty.

###  2.3 Bayesian hierarchical modelling and Max-and-Smooth

All our postprocessing methods are such that observations $y_{s,t}$ at spatial
locations $s=1,\dots,S$ and time instances $t=1,\dots,T$ are independently
distributed with some distribution $p(y_{s,t}|\boldsymbol{\theta}_{s})$ that
is conditional on spatially varying (but time-constant) postprocessing
parameters $\boldsymbol{\theta}_{s}=(\alpha_{s},\beta_{s},\dots)^{\prime}$.
The dependency of $p(y_{s,t}|\boldsymbol{\theta}_{s})$ on a vector of
covariates $\boldsymbol{z}_{s,t}$ is understood, and not written out
explicitly. Let observations at all locations and times be collected in the
vector $\boldsymbol{y}$ and all postprocessing parameters collected in the
vector $\boldsymbol{\theta}$. Conditional independence then implies

$\log p(\boldsymbol{y}|\boldsymbol{\theta})=\sum_{s=1}^{S}\sum_{t=1}^{T}\log
p(y_{s,t}|\boldsymbol{\theta}_{s}).$ (12)

For Bayesian inference of $\boldsymbol{\theta}$ we assume
$\boldsymbol{\theta}$ to have a prior normal distribution with mean vector
$\boldsymbol{\mu}_{\theta}$ and precision matrix $Q_{\theta}$, i.e., the prior
log-density of $\boldsymbol{\theta}$ can be written as

$\log
p(\boldsymbol{\theta})=const-\frac{1}{2}(\boldsymbol{\theta}-\boldsymbol{\mu}_{\theta})^{\prime}Q_{\theta}(\boldsymbol{\theta}-\boldsymbol{\mu}_{\theta}),$
(13)

where $const$ denotes a generic additive constant that does not depend on
$\boldsymbol{\theta}$. Defining

$f_{s}(\boldsymbol{\theta}_{s})=\sum_{t=1}^{T}\log
p(y_{s,t}|\boldsymbol{\theta}_{s}),$ (14)

the posterior log-density of $\boldsymbol{\theta}$ is given by

$\log
p(\boldsymbol{\theta}|\boldsymbol{y})=const+\sum_{s=1}^{S}f_{s}(\boldsymbol{\theta}_{s})-\frac{1}{2}(\boldsymbol{\theta}-\boldsymbol{\mu}_{\theta})^{\prime}Q_{\theta}(\boldsymbol{\theta}-\boldsymbol{\mu}_{\theta}).$
(15)

To make posterior inference of $\boldsymbol{\theta}$ analytically tractable,
we Taylor-approximate $f_{s}(\boldsymbol{\theta}_{s})$ to second order around
its mode $\hat{\boldsymbol{\theta}}_{s}$

$f_{s}(\boldsymbol{\theta}_{s})\approx\tilde{f}_{s}(\boldsymbol{\theta}_{s})=const-\frac{1}{2}(\boldsymbol{\theta}_{s}-\hat{\boldsymbol{\theta}}_{s})^{\prime}\hat{J}_{s}(\boldsymbol{\theta}-\hat{\boldsymbol{\theta}}_{s})$
(16)

where
$\hat{\boldsymbol{\theta}}_{s}=\mathrm{argmax}f_{s}(\boldsymbol{\theta}_{s})$
is the vector of maximum likelihood estimates (MLEs) at spatial location $s$,
and $\hat{J}_{s}$ is the observed information matrix at location $s$ evaluated
at the local MLEs $\hat{\boldsymbol{\theta}}_{s}$. Denote the combined vector
of MLEs at all spatial locations by $\hat{\boldsymbol{\theta}}$, and the
combined observed information matrix
$\hat{J}=J(\hat{\boldsymbol{\theta}})=-\frac{\partial^{2}}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{\prime}}\sum_{s=1}^{S}f_{s}(\boldsymbol{\theta}_{s})|_{\hat{\boldsymbol{\theta}}}$.
If the elements of the combined parameter vector $\boldsymbol{\theta}$ are
ordered as
$\boldsymbol{\theta}^{\prime}=(\boldsymbol{\theta}_{1}^{\prime},\dots,\boldsymbol{\theta}_{S}^{\prime})$,
then the combined information matrix is block-diagonal as
$\hat{J}=\mathrm{bdiag}(\hat{J}_{1},\dots,\hat{J}_{S})$. Under the second-
order Taylor approximation of $f_{s}(\boldsymbol{\theta}_{s})$ we have an
approximation of the posterior distribution of $\boldsymbol{\theta}$ given by

$\displaystyle\log p(\boldsymbol{\theta}|\boldsymbol{y})\approx$
$\displaystyle\log\tilde{p}(\boldsymbol{\theta}|\boldsymbol{y})$ (17)
$\displaystyle=const-\frac{1}{2}(\boldsymbol{\theta}-\hat{\boldsymbol{\theta}})^{\prime}\hat{J}(\boldsymbol{\theta}-\hat{\boldsymbol{\theta}})-\frac{1}{2}(\boldsymbol{\theta}-\boldsymbol{\mu}_{\theta})^{\prime}Q_{\theta}(\boldsymbol{\theta}-\boldsymbol{\mu}_{\theta})$
$\displaystyle=const-\frac{1}{2}\boldsymbol{\theta}^{\prime}(Q_{\theta}+\hat{J})\boldsymbol{\theta}+(Q_{\theta}\boldsymbol{\mu}_{\theta}+\hat{J}\hat{\boldsymbol{\theta}})^{\prime}\boldsymbol{\theta}$

which is the log-density of a multivariate normal distribution in canonical
form (Rue and Held, 2005), and hence the posterior distribution of
$\boldsymbol{\theta}$ is approximately multivariate normal with conditional
expectation

$E(\boldsymbol{\theta}|\boldsymbol{y})=(Q_{\theta}+\hat{J})^{-1}(Q_{\theta}\boldsymbol{\mu}_{\theta}+\hat{J}\hat{\boldsymbol{\theta}})$
(18)

and conditional variance matrix

$Var(\boldsymbol{\theta}|\boldsymbol{y})=(Q_{\theta}+\hat{J})^{-1}.$ (19)

The Gaussian approximation of the posterior is the central approximation
underlying the Integrated Nested Laplace Approximations method (INLA, Rue et
al., 2009) used by Möller et al. (2015) for spatial inference of
postprocessing parameters in a Bayesian hierarchical model. Their spatial
prior is a stochastic partial differential equation (SPDE) which corresponds
to a specific choice of the prior precision matrix $Q_{\theta}$. Our spatial
prior and corresponding choice of $Q_{\theta}$ is introduced in Section 2.4.

We now present an alternative derivation of Eq. 17 which makes the connection
between the approximate Bayesian hierarchical modelling approach with R-INLA
used by Möller et al. (2015), and the parameter smoothing approach for spatial
postprocessing proposed by Kharin et al. (2017). According to asymptotic
likelihood theory (Schervish, 2012) the MLE $\hat{\boldsymbol{\theta}}_{s}$
has an asymptotic normal (sampling) distribution centered on the “true”
parameter value $\boldsymbol{\theta}_{s}$ and with asymptotic precision given
by the expected information
$I_{s}=-E\left[\frac{\partial^{2}}{\partial\boldsymbol{\theta}_{s}\partial\boldsymbol{\theta}_{s}^{\prime}}\log
p(\boldsymbol{y}_{s}|\boldsymbol{\theta}_{s})\right]$. Since the observed
information $\hat{J}_{s}$ is often a more appropriate measure of the precision
of $\hat{\boldsymbol{\theta}}_{s}$ at estimating $\boldsymbol{\theta}_{s}$
(Efron and Hinkley, 1978) we may approximate the sampling distribution by

$\hat{\boldsymbol{\theta}}_{s}\overset{\cdot}{\sim}N(\boldsymbol{\theta}_{s},\hat{J}_{s}^{-1}).$
(20)

We interpret the MLEs $\hat{\boldsymbol{\theta}}_{s}$ as “noisy measurements”
of the “true” parameters $\boldsymbol{\theta}_{s}$, with normally distributed
measurement error with zero mean and variance $\hat{J}_{s}^{-1}$. Assuming
that the MLEs at different locations are independent we thus have

$\log
p(\hat{\boldsymbol{\theta}}|\boldsymbol{\theta})=const-\frac{1}{2}(\hat{\boldsymbol{\theta}}-\boldsymbol{\theta})^{\prime}\hat{J}(\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}).$
(21)

To infer $\boldsymbol{\theta}$ from the observed values of
$\hat{\boldsymbol{\theta}}$ we specify the same multivariate normal prior for
$\boldsymbol{\theta}$ as in Eq. 13. Then we have for the posterior of
$\boldsymbol{\theta}$

$\displaystyle\log p(\boldsymbol{\theta}|\hat{\boldsymbol{\theta}})$
$\displaystyle=const+\log
p(\hat{\boldsymbol{\theta}}|\boldsymbol{\theta})+\log p(\boldsymbol{\theta})$
(22)
$\displaystyle=const-\frac{1}{2}\boldsymbol{\theta}^{\prime}(Q_{\theta}+\hat{J})\boldsymbol{\theta}+(Q_{\theta}\boldsymbol{\mu}_{\theta}+\hat{J}\hat{\boldsymbol{\theta}})^{\prime}\boldsymbol{\theta}$
(23)

which implies

$\displaystyle E(\boldsymbol{\theta}|\hat{\boldsymbol{\theta}})$
$\displaystyle=(Q_{\theta}+\hat{J})^{-1}(Q_{\theta}\boldsymbol{\mu}_{\theta}+\hat{J}\hat{\boldsymbol{\theta}})\text{,
and}$ (24) $\displaystyle Var(\boldsymbol{\theta}|\hat{\boldsymbol{\theta}})$
$\displaystyle=(Q_{\theta}+\hat{J})^{-1}.$

The approximate posterior of $\boldsymbol{\theta}$ given
$\hat{\boldsymbol{\theta}}$ is identical to the approximate posterior of
$\boldsymbol{\theta}$ given $\boldsymbol{y}$ in Eq. 17. We have thus shown
that a normal approximation of the likelihood function in a Bayesian
hierarchical model with a multivariate normal prior results in the same
posterior as a model that assumes that the MLEs $\hat{\boldsymbol{\theta}}$
are “noisy measurements” of the “true” parameters $\boldsymbol{\theta}$, with
error (co-)variances derived from their asymptotic sampling distribution.
Since Eq. 20 models the MLEs as “truth plus error”, the inferred posterior
expectations of $\boldsymbol{\theta}$ will usually appear spatially smoother
than the MLEs $\hat{\boldsymbol{\theta}}$. The inference of the posterior
expectation of $\boldsymbol{\theta}$ can thus be regarded as a smoothing
applied to the local MLEs, which is also what Kharin et al. (2017) proposed to
use spatial information more effectively in postprocessing of gridded data.

In summary, to approximate spatial postprocessing by a Bayesian hierarchical
model we can proceed in a two-steps: First we postprocess forecasts locally,
by calculating the MLEs and the observed information matrix at each location
separately. Secondly we “post-postprocess” the MLEs by spatially smoothing
them, under a measurement error model using a spatially correlated prior and
measurement error (co-)variances derived from the observed information.
Hrafnkelsson et al. (2021) hence dubbed that approximation “Max-and-Smooth”.
While yielding identical results, our derivation of Max-and-Smooth differs
slightly from the original. Here, we motivate the measurement error model via
the asymptotic sampling distribution of the MLEs. Hrafnkelsson et al. (2021)
simply “reverse-engineered” the measurement error model from the Laplace
approximation of the posterior, and showed that it yields identical results.

###  2.4 The RW2D spatial prior

The 2-dimensional random walk (RW2D, Rue and Held, 2005) is a first order
intrinsic autoregressive model on a lattice indexed by row $i=1,\dots,N_{i}$
and column $j=1,\dots,N_{j}$. An RW2D random field
$W\in\mathbb{R}^{N_{i}\times N_{j}}$ is defined recursively by assuming an
independent normal distribution of 2-dimensional increments,

$W_{i,j}-\frac{1}{4}\left(W_{i-1,j}+W_{i+1,j}+W_{i,j-1}+W_{i,j+1}\right)\sim
N\left(0,\frac{1}{\kappa}\right).$ (25)

The relationship between $W_{i,j}$ and the average over its nearest neighbors
induces spatial correlation. The higher the precision parameter $\kappa$ the
closer $W_{i,j}$ is (on average) to its neighbors, and hence $\kappa$ controls
the smoothness of the random field $W$. For more background on the RW2D see
Besag and Mondal (2005), who discuss the relationship between intrinsic
autoregressions and the de Wijs process, and their examples include the RW2D
as a special case.

Let $W\in\mathbb{R}^{N_{i}\times N_{j}}$ be a realisation of an RW2D on a
$N_{i}\times N_{j}$ grid, and $E\in\mathbb{R}^{N_{i}\times N_{j}}$ a matrix of
independent standard normally distributed random variates. Then the RW2D can
be written as the Sylvester equation

$D_{N_{i}}W+WD_{N_{j}}=\kappa^{-1/2}E,$ (26)

where $D_{n}$ denotes the tri-diagonal $n\times n$ matrix

$D_{n}=\frac{1}{4}\begin{pmatrix}1&-1\\\ -1&2&-1\\\ &\ddots&\ddots&\ddots\\\
&&-1&2&-1\\\ &&&-1&1\end{pmatrix}.$ (27)

Let the $N_{i}N_{j}$-vectors $\boldsymbol{w}$ and $\boldsymbol{e}$ be the
vectorisations of $W$ and $E$ obtained by column-stacking. Then Eq. 26 can be
written as

$\kappa^{1/2}D\boldsymbol{w}=\boldsymbol{e}$ (28)

with

$D=D_{N_{i}}\otimes I_{N_{j}}+I_{N_{i}}\otimes D_{N_{j}},$ (29)

where $I_{n}$ denotes the $n\times n$ identity matrix (Shores, 2007). It
follows that $\boldsymbol{w}$ is a normally distributed random vector with
mean zero and precision matrix

$Q_{RW2D}=\kappa R_{RW2D}=\kappa D^{\prime}D$ (30)

(Rue and Held, 2005). We will use $Q_{RW2D}$ as the spatial prior precision
$Q_{\theta}$ to infer postprocessing parameters throughout this paper.

The $S\times S$ matrix $Q_{RW2D}$ (where $S=N_{i}N_{j}$) is rank-deficient
with $\mathrm{rank}(Q_{RW2D})=S-1$, and hence $Q_{RW2D}$ does not have a well-
defined inverse. Marginal variances and covariances of the RW2D process are
thus ill-defined, and numerical sampling from the RW2D process problematic. In
our application, $Q_{RW2D}$ is used merely as a prior to infer the posterior
mean and variance of $\boldsymbol{\theta}$ (Eq. 18 and Eq. 19). These
posterior quantities depend on inverses of $(Q_{RW2D}+\hat{J})$ which is
usually of full rank, and hence the rank deficiency of $Q_{RW2D}$ is benign.

For joint spatial inference of multiple postprocessing parameters (e.g.
$\alpha^{LR}$ and $\beta^{LR}$ in logistic regression), we specify independent
spatial RW2D priors for them. If the elements of the full parameter vector
$\boldsymbol{\theta}$ are ordered as
$\boldsymbol{\theta}=(\alpha_{1},\dots,\alpha_{S},\beta_{1},\dots,\beta_{S},\dots)^{\prime}$,
and each of $\boldsymbol{\alpha},\boldsymbol{\beta},\dots$ have independent
RW2D priors, then the joint prior for $\boldsymbol{\theta}$ is a multivariate
normal distribution with mean zero and block-diagonal precision matrix

$Q_{\theta}=\begin{pmatrix}Q_{\alpha}&&\\\ &Q_{\beta}&\\\
&&\ddots\end{pmatrix}=\begin{pmatrix}\kappa_{\alpha}R_{Rw2D}&&\\\
&\kappa_{\beta}R_{RW2D}&\\\ &&\ddots\end{pmatrix}.$ (31)

The precision parameters
$\boldsymbol{\kappa}=(\kappa_{\alpha},\kappa_{\beta},\dots)^{\prime}$ are
hyperparameters.

The matrices $Q_{RW2D}$ and $\hat{J}$ are sparse matrices. Sparsity allows for
efficient numerical solutions of linear systems such as Eq. 18, and efficient
random sampling which can be used to calculate approximate marginal variances
from Eq. 19 (Rue and Held, 2005). In this paper, we use sparse matrix methods
implemented in the R package Matrix (Bates et al., 2022).

###  2.5 Inferring hyperparameters

For spatial inference of postprocessing parameters with a fixed prior
precision matrix $Q_{\theta}$ we have to specify the precision parameter
$\kappa$ for each parameter, which controls the amount of smoothing when
inferring $\boldsymbol{\theta}$ from $\hat{\boldsymbol{\theta}}$. Let the
vectors $\boldsymbol{\kappa}$ and $\boldsymbol{\theta}$ denote vectors of,
respectively, hyperparameters and postprocessing parameters. We momentarily
write the spatial prior for $\boldsymbol{\theta}$ as conditional on
$\boldsymbol{\kappa}$, i.e. $p(\boldsymbol{\theta}|\boldsymbol{\kappa})$ and
specify a prior distribution $p(\boldsymbol{\kappa})$. The posterior log-
density of the hyperparameters $\boldsymbol{\kappa}$, given the MLEs
$\hat{\boldsymbol{\theta}}$, is given by

$\log p(\boldsymbol{\kappa}|\hat{\boldsymbol{\theta}})=\log
p(\hat{\boldsymbol{\theta}}|\boldsymbol{\theta})+\log
p(\boldsymbol{\theta}|\boldsymbol{\kappa})+\log p(\boldsymbol{\kappa})-\log
p(\boldsymbol{\theta}|\hat{\boldsymbol{\theta}},\boldsymbol{\kappa})-\log
p(\hat{\boldsymbol{\theta}})$ (32)

where $p(\hat{\boldsymbol{\theta}}|\boldsymbol{\theta})$ follows from the
normally distributed measurement error model of the MLEs, and the posterior
$p(\boldsymbol{\theta}|\hat{\boldsymbol{\theta}},\boldsymbol{\kappa})$ is
approximated by a normal distribution as in Eq. 24. The term
$p(\hat{\boldsymbol{\theta}})$ is a constant that does not depend on
$\boldsymbol{\kappa}$, and the value of $\boldsymbol{\theta}$ when evaluating
the rhs of Eq. 32 can be set arbitrarily.

For the analyses in this paper we infer $\boldsymbol{\theta}$ by empirical
Bayes, i.e. setting the precision hyperparameters
$\boldsymbol{\kappa}=(\kappa_{\alpha},\kappa_{\beta},\dots)^{\prime}$ to the
fixed values that maximise Eq. 32 with respect to $\boldsymbol{\kappa}$.
Ignoring any terms that do not depend on $\boldsymbol{\kappa}$, we thus
estimate $\boldsymbol{\kappa}$ by the point estimator

$\displaystyle\hat{\boldsymbol{\kappa}}$
$\displaystyle=\underset{\boldsymbol{\kappa}}{\mathrm{argmax}}\Big{[}-\frac{1}{2}(\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*})^{\prime}\hat{J}(\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}^{*})+\frac{1}{2}\log\det
Q_{\theta}$ (33)
$\displaystyle\quad\qquad\qquad-\frac{1}{2}{\boldsymbol{\theta}^{*}}^{\prime}Q_{\theta}\boldsymbol{\theta}^{*}+\log
p(\boldsymbol{\kappa})-\frac{1}{2}\log\det(Q_{\theta}+\hat{J})\Big{]}.$

where
$\boldsymbol{\theta}^{*}=E(\boldsymbol{\theta}|\hat{\boldsymbol{\theta}},\boldsymbol{\kappa})$
given by Eq. 24. Since $Q_{\theta}$ is block-diagonal (Eq. 31) and $R_{RW2D}$
has rank $S-1$ and is independent of $\boldsymbol{\kappa}$, the log-
determinant of $Q_{\theta}$ is given by $\log\det
Q_{\theta}=const+(S-1)(\log\kappa_{\alpha}+\log\kappa_{\beta}+\dots)$ (Rue and
Held, 2005). We use independent Exponential prior distributions with rate
$5\times 10^{-5}$ for the elements of $\boldsymbol{\kappa}$, hence

$\log p(\boldsymbol{\kappa})=const-5\times
10^{-5}(\kappa_{\alpha}+\kappa_{\beta}+\dots).$ (34)

To find $\hat{\boldsymbol{\kappa}}$ we use numerical optimisation in the R
function optim with default settings.

The obtained value of $\hat{\boldsymbol{\kappa}}$ is plugged into $Q_{\theta}$
to calculate the smoothed postprocessing via Eq. 18. Our empirical Bayes
approach ignores the estimation uncertainty due to inferring
$\boldsymbol{\kappa}$. To assess the effect of this simplification we compare
our results to parameter estimates obtained with R-INLA (cf. Sec. 2.6) which
approximates a fully Bayesian inference where $\boldsymbol{\kappa}$ is
integrated out rather than fixed.

###  2.6 R-INLA

The use of multivariate normal priors with sparse precision matrices, and the
normal approximation of the likelihood we used in Eq. 16 are also the main
ingredients used in the R-INLA package (Rue et al., 2022) for large scale
Bayesian inference. R-INLA has the option of inferring the posterior of the
latent field $\boldsymbol{\theta}$ by empirical Bayes, i.e.
$p(\boldsymbol{\theta}|\boldsymbol{y})=p(\boldsymbol{\theta}|\boldsymbol{y},\hat{\boldsymbol{\kappa}})$
or by numerically integrating over the posterior of hyperparameters, i.e.,
solving $p(\boldsymbol{\theta}|\boldsymbol{y})=\int
p(\boldsymbol{\theta}|\boldsymbol{y},\boldsymbol{\kappa})p(\boldsymbol{\kappa}|\boldsymbol{y})d\boldsymbol{\kappa}$.
There are additional refinements to the various approximations as detailed in
Rue et al. (2009), but essentially our approach is very similar to R-INLA. A
review of R-INLA’s spatial modelling capabilities can be found in Bakka et al.
(2018), and coding tutorials and examples are available in
https://www.r-inla.org/examples-tutorials222Last accessed 6 July 2022.,
Blangiardo and Cameletti (2015) and Gómez-Rubio (2020).

Möller et al. (2015) use R-INLA to infer
$p(\boldsymbol{\theta}|\boldsymbol{y})$, by specifying a Bayesian hierarchical
model for the meteorological observations $\boldsymbol{y}$ directly, without
the intermediate step of calculating MLEs first. While it is possible to
specify RW2D and many other spatial priors in R-INLA, it cannot define all
likelihood functions commonly used for forecast postprocessing. For example,
while Bayesian linear regression is readily implemented in R-INLA, it is not
possible333as of version 22.01.12 to specify a spatial prior for the variance
parameter, or make the variance dependent on covariates as in NGR. R-INLA is
thus not applicable in many spatial post-processing applications.

Bayesian Logistic Regression with spatial priors is, however, implemented in
R-INLA, so we can use R-INLA as a reference implementation. We can also assess
the effects of our empirical Bayes approach of fixing hyperparameters at
$\hat{\boldsymbol{\kappa}}$ compared to integrating them out. R-INLA includes
corrections to improve the normal approximation of the likelihood function,
whose benefits we will also be able to assess in a comparison of spatial
logistic regression with R-INLA.

Under an additional approximation of the observation matrix, R-INLA can be
used to perform the smoothing step of Max-and-Smooth. Consider the diagonal
approximation of $\tilde{J}$ for $\hat{J}$

$\tilde{J}=\left[\mathrm{diag}(\hat{J}^{-1})\right]^{-1},$ (35)

i.e. the variance matrix of the MLEs is a diagonal matrix with elements equal
to their marginal asymptotic variances and covariances between MLEs are
ignored. In the case of MOS that approximation is exact because the MLEs are
asymptotically independent. When $Q_{\theta}$ is block-diagonal, with one
block per parameter, and the information matrix is diagonal (either by
construction as in MOS, or by approximation), then the posterior mean
$(Q_{\theta}+\tilde{J})^{-1}\tilde{J}\hat{\boldsymbol{\theta}}$ can be
decomposed into updates of individual parameters, i.e.
$E(\boldsymbol{\alpha}|\hat{\boldsymbol{\alpha}})=(Q_{\alpha}+\tilde{J}_{\alpha})^{-1}\tilde{J}_{\alpha}\hat{\boldsymbol{\alpha}}$,
etc, where $Q_{\alpha}$ is the RW2D prior precision for $\boldsymbol{\alpha}$.
The diagonal approximation allows smoothing the MLEs individually and
independently, which is straightforward with R-INLA. The approximation Eq. 35
was used in Jóhannesson et al. (2022) for model selection, and Lovegrove and
Siegert (2022) used it to apply INLA for spatial postprocessing with MOS.
Example R code for smoothing a spatial field of MLEs with known sampling
variances with R-INLA is provided in the Supplementary Material.

###  2.7 Model evaluation

We illustrate and discuss the performance of the proposed methods by
postprocessing the forecast data presented in Section 2.1. In all analyses we
first apply local postprocessing by maximum likelihood estimation separately
at each grid point and lead time. We apply spatial smoothing to the
2-dimensional fields of local MLEs using Max-and-Smooth with independent RW2D
priors. The quality of postprocessing models is evaluated and compared by
deterministic and probabilistic verification metrics.

We are working with complete data on a rectangular grid throughout. There are
a total of $S$ grid points, indexed by the single spatial index $s=1,\dots,S$.
At each grid point $s$ and lead time $\tau$, a number of $T$ forecast-
observation pairs is available at verification times indexed by $t=1,\dots,T$.
Depending on context, we average verification metrics either over both,
spatial coordinate $s$ and verification time $t$, or, to highlight spatial
differences of performance, only over verification time $t$. Verification
scores are always calculated separately at each lead time $\tau$, and hence
the dependency on lead time $\tau$ is suppressed in mathematical notation.

The accuracy of postprocessed temperature forecasts (i.e. , expectations of
the postprocessed predictive distributions) are evaluated by calculating the
mean squared errors (MSE) (Déqué, 2012). Let $\mu_{s,t}$ denote the mean of
the forecast distribution at spatial location $s$ and verification time $t$,
and let $y_{s,t}$ be the corresponding verifying observation. The MSE is the
space-time averaged squared prediction error

$MSE=\frac{1}{ST}\sum_{s=1}^{S}\sum_{t=1}^{T}(y_{s,t}-\mu_{s,t})^{2}.$ (36)

Forecast probabilities $p_{s,t}$ of precipitation exceedances (binary yes/no,
encoded as $y_{s,t}=1$ or $0$, respectively) are evaluated by the Brier score
(Brier, 1950). The average Brier score is given by the mean squared difference
between forecast probabilities $p_{s,t}=P(y_{s,t}=1)$ and binary event
indicators

$BS=\frac{1}{ST}\sum_{s=1}^{S}\sum_{t=1}^{T}(y_{s,t}-p_{s,t})^{2}.$ (37)

The Brier score is a proper scoring rule that rewards truthful representations
of the probability distribution of the target event frequency Broecker (2012).
Like the MSE it is negatively oriented, assigning lower values to better
forecasts.

Probabilistic temperature forecast, issued as full probability density
functions (pdfs) over all possible outcomes are evaluated by the logarithmic
score (also Logscore or Ignorance, Bernardo, 1979; Roulston and Smith, 2002).
Let $f_{s,t}(x)$ denote the forecast pdf at grid point $s$ and verification
time $t$, and $y_{s,t}$ the real-valued verifying observation. The Logscore is
given by

$LS=\frac{1}{ST}\sum_{s=1}^{S}\sum_{t=1}^{T}\left[-\log
f_{s,t}(y_{s,t})\right].$ (38)

The Logscore is proper and negatively oriented.

Our goal is to assess whether, and how much, parameter smoothing can improve
forecasts of future observations in operational settings. To assess out-of-
sample performance, and avoid results that are biased by overfitting, we apply
leave-one-out cross validation (Hastie et al., 2009). Specifically, we leave
each individual verification time $t=1,\dots,T$ out in turn, estimate
parameters (unsmoothed and smoothed) and hyperparameters on the data set
comprising the remaining $T-1$ verification times, and calculate evaluation
metrics by predicting observations at the left-out verification time $t$. We
repeat this process $T$ times and average evaluation metrics over all $T$. In
our case, the relatively small number of verification times ($T=20$), as well
as the efficiency of the spatial smoothing method makes this “brute-force”
leave-one-out cross validation practically feasible, without having to resort
to approximate model selection critiera (Hastie et al., 2009; Piironen and
Vehtari, 2017).

##  3 Results

In this section we apply Max-and-Smooth to spatial postprocessing of
precipitation forecasts with Logistic Regression, and spatial postprocessing
of temperature forecasts with Model Output Statistics (MOS) and Nonhomogeneous
Gaussian Regression (NGR). The focus of the analyses is on the comparison of
out-of-sample accuracy and probabilistic skill of forecasts before and after
spatial smoothing of post-processing parameters.

###  3.1 Logistic regression

We fit Logistic Regression to estimate the probability of threshold exceedance
based on the total precipitation ensemble mean. We estimate postprocessing
parameters by maximum likelihood estimation as described in Sec. 2.2. We then
calculate posterior means of postprocessing parameters via Max-and-Smooth, by
applying Eq. 24, with hyperparameters estimated using Eq. 33. We also use the
R-INLA package to fit postprocessing parameters in a Bayesian logistic
regression model with spatial RW2D priors on the intercept and slope
parameters. The main difference between R-INLA and our Max-and-Smooth
implementation is that hyperparameters are integrated out, rather than plugged
in. We also evaluate the local climatological forecasts (the unconditional
mean rate of threshold exceedances) as a benchmark.

Figure 1: Parameter estimates of a logistic regression of the probability of
exceeding precipitation of $2.5$mm at 1 day lead time. Upper panels: Local
maximum-likelihood estimates of regression parameters. Lower Panel: Posterior
means obtained by Max-and-Smooth.

Figure 1 compares local MLEs with posterior means of postprocessing parameters
calculated with Max-and-Smooth. The MLEs fluctuate a lot spatially, and so
Bayesian inference with a spatially correlated prior has a strong smoothing
effect on the parameter estimates. After applying the smoothing step, only
large scale patterns of variability remain. A noteable result is that after
spatial smoothing all slope parameters of the regression are strictly
positive.

Figure 2: Brier scores (mean $\pm$ two standard deviations) of climatological
forecast, logistic regression using local MLEs, INLA, and Max-and-Smooth.

Figure 2 shows the spatially averaged leave-one-out Brier scores obtained by
local postprocessing, and obtained by spatial post-processing using Max-and-
Smooth and R-INLA. We note from Figure 2 that predicting the threshold
exceedance based on the ensemble mean seems to be a difficult forecasting
problem with generally low skill. This is indicated by the poor performance of
the local MLE logistic regression compared to climatology, even at short lead
times. However, forecast skill can be improved substantially by spatial
postprocessing. The absolute Brier Score improvement of Max-and-Smooth over
MLE is between $0.14$ and $0.18$, which corresponds to relative improvements
between $8.5\%$ and $11.1\%$. We further note that Max-and-Smooth achieves
similar forecast skill as INLA, which suggests that integrating out
hyperparameters has little effect compared to fixing $\boldsymbol{\kappa}$ at
a point estimate.

###  3.2 MOS

In this section we study the skill of temperature forecasts postprocessed with
MOS, and compare the performance of postprocessing parameters estimated with
MLE and Max-and-Smooth.

Figure 3: Spatial maps of MOS coefficients before and after smoothing. The
forecast lead time is 1 day.

Figure 3 show MOS parameter estimates obtained by MLE and Max-and-Smooth. In
all parameters most of the spatial variations (small and large scale) are
preserved after smoothing. Hence the postprocessed forecasts obtained by MLE
and Max-and-Smooth will not differ much, and we can only expect small
differences in forecast skill, if any. The slope parameters undergo more
smoothing than the intercept and log-variance parameters.

Figure 4: Scatter plots of (50 randomly sampled) Max-and-Smooth parameter
estimates over the corresponding MLEs. The half-widths of the horizontal lines
are equal to two standard deviations of the asymptotic sampling distribution
of the MLEs.

The scatter plots of MOS coefficients in Figure 4 explain why the slopes are
smoothed more than the intercepts and log-variances. The asymptotic errors of
the intercept MLEs are small, and hence the inferred value of $\alpha_{s}$ at
each location is highly constrained to be close to $\hat{\alpha}_{s}$,
resulting in very little variation during smoothing. On the other hand, the
asymptotic errors of the slopes are relatively large. Hence $\hat{\beta}_{s}$
is less informative about $\beta_{s}$, and the spatial model has more freedom
to move the parameters around, resulting in more smoothing.

Figure 5: Mean squared prediction error difference ($\pm$ 2 standard
deviations) between postprocessed forecasts using local MLE vs. using Max-and-
Smooth. Spatial parameter smoothing with Max-and-Smooth improves the MSE at
all lead times.

Figure 5 shows MSE differences between local and spatially postprocessed
temperature forecasts. Consistent with results in section 3.1, there is a
small, but systematic improvement in MSE at all lead times.

Figure 6: Improvement of mean squared prediction error of Max-and-Smooth vs.
MLE.

Figure 6 shows the spatial distribution of MSE differences for lead times 2,
5, and 8 days. At most locations the MSE differences are positive indicating a
consistent improvement of forecast accuracy due to Max-and-Smooth. Bigger
improvements in MSE are seen over land than over sea. There are a few
locations where MSE deteriorates when using Max-and-Smooth. But the typical
magnitudes of positive differences are greater than those of negative
differences, resulting in an overall average improvement.

Figure 7: Spatially averaged Logscore difference of forecasts postprocessed
with MOS, using MLE parameters and Max-and-Smooth parameters. Positive values
indicate improvements due to Max-and-Smooth.

Analysis of probabilistic forecast skill by Logscores shown in Figure 7
confirms the improvement seen in terms of MSE. The forecast distributions
obtained with Max-and-Smooth have better Logscores than forecasts obtained
with the unsmoothed MLEs. The improvement is between $0.004$ and $0.02$,
indicating that the Max-and-Smooth distributions assign between
$\exp(0.004)\approx 1.004$ and $\exp(0.02)\approx 1.02$ times more density to
observations than the MLE distributions.

###  3.3 Non-homogeneous Gaussian regression

Here use instantaneous ensemble means and variances to postprocess temperature
forecasts by NGR, as described in section 2.2. By letting the forecast
distribution be informed by the ensemble variance through the parameters
$\gamma_{s}$ and $\delta_{s}$, we allow narrower (wider) ensembles to generate
a narrower (wider) postprocessed forecast distribution. As before, we first
estimate grid-point wise NGR parameters by MLE, and subsequently smooth them
using Eq. 24, with smoothing parameters $\boldsymbol{\kappa}$ estimated with
Eq. 33. We also smooth postprocessing parameters individually with R-INLA,
using the diagonal approximation of the information matrix discussed in Sec.
2.6. Smoothing the MLEs independently is only an approximation to “proper”
Max-and-Smooth because covariances between MLEs are ignored.

Figure 8: Spatially averaged Logscores of NGR models using local MLEs, spatial
parameter smoothing with R-INLA assuming independence, and smoothing with Max-
and-Smooth accounting for parameter dependence. Smaller values indicate better
forecasts.

Figure 8 shows average Logscores as function of lead time of local MLE
forecasts, forecasts based on individually smoothed postprocessing parameters
with R-INLA, and forecasts based on jointly smoothed parameters using Max-and-
Smooth. Logscores generally increase with increasing forecast lead time, as
would be expected. Max-and-Smooth consistently outperforms both the local MLE
forecasts, and the independently smoothed INLA forecasts, highlighting the
benefit of spatial smoothing as well as the importance of the covariances
between MLEs. Interestingly, at lead times above 5 days, the independently
smoothed parameters perform worse than the unsmoothed MLEs. The improvements
in Logscores are between $0.09$ (lead time 10 days) and $0.15$ (lead time 1
day), indicating that applying Max-and-Smooth leads to forecast densities that
assign between $e^{0.09}\approx 1.09$ and $e^{0.15}\approx 0.16$ more density
to the observations than the local MLEs.

We further analysed probabilistic skill by the Continuous Ranked Probability
Score (CRPS, Matheson and Winkler, 1976), and forecast calibration by
Probability Integral Transform (PIT) histograms (Gneiting et al., 2007).
Details of these analyses can be seen in the Supplementary Material. In
summary, average CRPS of forecasts obtained with independently smoothed
parameters is consistently worse than with unsmoothed parameters, and
consistently better when parameters are postprocessed with Max-and-Smooth,
accounting for depdendency between parameters. PIT histograms obtained with
unsmoothed parameters appear U-shaped, but close to uniform after parameter
smoothing with either method, indicating considerable improvement of forecast
calibration.

Figure 9: Estimates of NGR variance parameters $\gamma$ and $\delta$, lead
time 1 day, obtained by MLE, smoothing with R-INLA assuming independence, and
Max-and-Smooth taking into account dependence between postprocessing
parameters.

Figure 9 shows the NGR variance parameters $\gamma$ and $\delta$ for the 1 day
lead time temperature forecasts, estimated by the three methods. The MLEs
appear spatially very noisy, and the effect of spatial smoothing by both
methods is large. The smoothed estimates obtained with R-INLA and obtained
with Max-and-Smooth look generally similar, but the Max-and-Smooth estimates
retain more spatial structure than the independent R-INLA estimates. Under the
diagonal approximation of the information matrix, a parameter at grid point
$i,j$ is adjusted only based on values of the same parameter at neighboring
grid points. In “full” Max-and-Smooth, which takes correlation between
parameters into account, the smoothed value of, say, $\delta$ at grid point
$i,j$ depends not only on values of $\delta$ at neighboring grid points, but
also on nearby values of $\alpha$, $\beta$ and $\gamma$. That dependency on
more information might explain the more detailed spatial structure of Max-and-
Smooth estimates compared to the independent R-INLA estimates. The effect is
large enough to lead to the notable improvements of forecast skill seen in
Fig. 8. We also note, the smoothed versions of the constant additive variance
offset, $e^{\gamma_{s}}$, are systematically greater than 1 over land, and
smaller than 1 over sea, indicating that ensemble variances exhibit more
underdispersion over land than over sea. That pattern is consistent across all
lead times (not shown).

##  4 Discussion

Spatial smoothing of postprocessing parameters has been previously proposed to
use spatial information more effectively (Kharin et al., 2017). We have shown
that parameter smoothing can be motivated by approximating Bayesian inference
of postprocessing parameters in a Bayesian hierarchical model with a spatially
correlated prior. Unlike Kharin et al. (2017), our approach explicitly
accounts for the estimation uncertainty of the local postprocessing
parameters. MLEs with smaller estimation uncertainty are smoothed less than
more uncertain parameters, and so considering estimation uncertainty is useful
to constrain the smoothing hyperparameters.

Max-and-Smooth applied to postprocessing leads to consistent improvements of
postprocessed forecasts in applications of logistic regression, linear
regression (MOS), and nonhomogenous Gaussian regression. Improvements were
demonstrated in terms of accuracy of the forecast mean, as measured by their
mean squared errors, probabilistic skill as measured by Brier scores,
Logscores, and CRPS, and probabilistic calibration as measured by PIT
histograms. Overall, we deem the evidence that spatial postprocessing is
beneficial convincing, and did therefore not carry out any significance
testing.

The proposed approach is simple to apply, as it is based on a straightforward
linear transformation of local parameter estimates. As such it operates on
results that practitioners generate routinely. Spatial postprocessing can be
regarded as a “post-postprocessing” step, hence maintaining a breakdown of
postprocessing steps into local postprocessing by calculating MLEs (“Max” …)
and subsequent spatial postprocessing of the local MLEs (… “and Smooth”). We
hope that such a clear separation of methodologies will allow for easier
adoption of the method in operational settings.

There are several directions for further study and application of the proposed
methodology. Max-and-Smooth involves approximations and simplifications whose
impact on reliability, accuracy, and skill of postprocessed forecasts can be
studied further. We have not implemented a fully Bayesian framework in which
parametric uncertainty in hyperparameters and postprocessing parameters is
propagated into the predictive distributions. Such an extension is
straightforward and might yield further improvements, although the comparison
between INLA and Max-and-Smooth in Sec. 3.1 suggest that these effects are
likely small. Max-and-Smooth presented here is a spatial postprocessing method
at the parameter level, and as such does not produce forecasts that are
spatially consistent with observations. Möller et al. (2015) showed that an
additional postprocessing step can be applied to independent samples from
postprocessed forecast distributions to recover spatial consistency, and such
an approach would also be straightforward to apply in our framework. For a
preliminary application of combining an empirical copula approach and Max-and-
Smooth to generate spatially coherent Generalised Extreme Value distributions,
see Jóhannesson et al. (2022). Gneiting et al. (2005) estimated NGR parameters
by minimum-CRPS estimation, rather than MLE. The theory to derive Max-and-
Smooth as an approximation of inference in a Bayesian hierarchical model
presented here requires use of the MLEs. But on a purely pragmatic level,
spatial smoothing of minimum-CRPS estimators is certainly possible, especially
considering the asymptotic Normality of such estimators (Yuen and Stoev,
2014). Lastly, spatial modelling allows for estimation of postprocessing
parameters in locations where no observation data are available. A spatial
modelling approach could thus be used to avoid the use of reanalysis data and
apply postprocessing directly to “real observations”, using (spatially sparse)
station data.

## Acknowledgments

In addition to the software packages cited in the article, data analyses and
visualisations relied on the R packages tidyverse (Wickham et al., 2019) and
Rnaturalearth (South, 2017). We thank Tommy Irons and Annette Möller for
helpful discussions and encouraging feedback.

## References

  * Bakka et al. [2018] H. Bakka, H. Rue, G.-A. Fuglstad, A. Riebler, D. Bolin, J. Illian, E. Krainski, D. Simpson, and F. Lindgren.  Spatial modeling with R-INLA: A review.  _Wiley Interdisciplinary Reviews: Computational Statistics_ , 10(6):e1443, 2018.
  * Bates et al. [2022] D. Bates, M. Maechler, and M. Jagan.  _Matrix: Sparse and Dense Matrix Classes and Methods_ , 2022.  URL https://CRAN.R-project.org/package=Matrix.  R package version 1.4-1.
  * Bernardo [1979] J. M. Bernardo.  Expected information as expected utility.  _The Annals of Statistics_ , pages 686–690, 1979.
  * Berrocal et al. [2008] V. J. Berrocal, A. E. Raftery, and T. Gneiting.  Probabilistic quantitative precipitation field forecasting using a two-stage spatial model.  _The Annals of Applied Statistics_ , 2(4):1170–1193, 2008.
  * Besag and Mondal [2005] J. Besag and D. Mondal.  First-order intrinsic autoregressions and the de Wijs process.  _Biometrika_ , 92(4):909–920, 2005.
  * Blangiardo and Cameletti [2015] M. Blangiardo and M. Cameletti.  _Spatial and spatio-temporal Bayesian models with R-INLA_.  John Wiley & Sons, 2015.
  * Brier [1950] G. W. Brier.  Verification of forecasts expressed in terms of probability.  _Monthly weather review_ , 78(1):1–3, 1950.
  * Broecker [2012] J. Broecker.  Probability forecasts.  In I. T. Jolliffe and D. B. Stephenson, editors, _Forecast Verification: A Practitioner’s Guide in Atmospheric Science_ , pages 118–139. John Wiley & Sons, Ltd., 2012.
  * Clark et al. [2004] M. Clark, S. Gangopadhyay, L. Hay, B. Rajagopalan, and R. Wilby.  The Schaake shuffle: A method for reconstructing space–time variability in forecasted precipitation and temperature fields.  _Journal of Hydrometeorology_ , 5(1):243–262, 2004.
  * Déqué [2012] M. Déqué.  Deterministic forecasts of continuous variables.  In I. T. Jolliffe and D. B. Stephenson, editors, _Forecast Verification: A Practitioner’s Guide in Atmospheric Science_ , pages 77–94. John Wiley & Sons, Ltd., 2012.
  * Efron and Hinkley [1978] B. Efron and D. V. Hinkley.  Assessing the accuracy of the maximum likelihood estimator: Observed versus expected fisher information.  _Biometrika_ , 65(3):457–483, 1978.
  * Glahn and Lowry [1972] H. R. Glahn and D. A. Lowry.  The use of model output statistics (MOS) in objective weather forecasting.  _Journal of Applied Meteorology_ , 11(8):1203–1211, 1972.
  * Gneiting et al. [2005] T. Gneiting, A. E. Raftery, A. H. Westveld, and T. Goldman.  Calibrated probabilistic forecasting using ensemble model output statistics and minimum crps estimation.  _Mon. Wea. Rev._ , 133(5):1098–1118, 2005.
  * Gneiting et al. [2007] T. Gneiting, F. Balabdaoui, and A. E. Raftery.  Probabilistic forecasts, calibration and sharpness.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 69(2):243–268, 2007.
  * Gómez-Rubio [2020] V. Gómez-Rubio.  _Bayesian inference with INLA_.  CRC Press, 2020.  Available online at https://becarioprecario.bitbucket.io/inla-gitbook/ (Last accessed on 6 July 2022).
  * Hastie et al. [2009] T. Hastie, R. Tibshirani, J. H. Friedman, and J. H. Friedman.  _The Elements of Statistical Learning: Data Mining, Inference, and Prediction_ , volume 2.  Springer, 2009.  Online: hastie.su.domains/ElemStatLearn.
  * Hersbach et al. [2020] H. Hersbach, B. Bell, P. Berrisford, S. Hirahara, A. Horányi, J. Muñoz-Sabater, J. Nicolas, C. Peubey, R. Radu, D. Schepers, et al.  The ERA5 global reanalysis.  _Quarterly Journal of the Royal Meteorological Society_ , 146(730):1999–2049, 2020.
  * Hrafnkelsson et al. [2021] B. Hrafnkelsson, S. Siegert, R. Huser, H. Bakka, and Á. V. Jóhannesson.  Max-and-Smooth: a two-step approach for approximate Bayesian inference in latent Gaussian models.  _Bayesian Analysis_ , 16(2):611–638, 2021.
  * Jóhannesson et al. [2022] Á. V. Jóhannesson, S. Siegert, R. Huser, H. Bakka, and B. Hrafnkelsson.  Approximate bayesian inference for analysis of spatiotemporal flood frequency data.  _The Annals of Applied Statistics_ , 16(2):905–935, 2022.
  * Kharin et al. [2017] V. Kharin, W. Merryfield, G. Boer, and W.-S. Lee.  A postprocessing method for seasonal forecasts using temporally and spatially smoothed statistics.  _Monthly Weather Review_ , 145(9):3545–3561, 2017\.  doi: https://doi.org/10.1175/MWR-D-16-0337.1.
  * Lovegrove and Siegert [2022] J. Lovegrove and S. Siegert.  Improving numerical weather forecasts by Bayesian hierarchical modelling.  In B. Hrafnkelsson, editor, _Statistical Modeling Using Bayesian Latent Gaussian models - with applications in geophysics and environmental sciences (to appear)_. Springer, 2022.
  * Matheson and Winkler [1976] J. E. Matheson and R. L. Winkler.  Scoring rules for continuous probability distributions.  _Management Science_ , 22(10):1087–1096, 1976\.
  * Möller et al. [2015] A. Möller, T. L. Thorarinsdottir, A. Lenkoski, and T. Gneiting.  Spatially adaptive, bayesian estimation for probabilistic temperature forecasts.  _arXiv preprint arXiv:1507.05066_ , 2015.
  * Piironen and Vehtari [2017] J. Piironen and A. Vehtari.  Comparison of bayesian predictive methods for model selection.  _Statistics and Computing_ , 27(3):711–735, 2017\.
  * Roulston and Smith [2002] M. S. Roulston and L. A. Smith.  Evaluating probabilistic forecasts using information theory.  _Monthly Weather Review_ , 130(6):1653 – 1660, 2002.
  * Rue and Held [2005] H. Rue and L. Held.  _Gaussian Markov Random Fields: Theory and Applications_.  CRC press, 2005.
  * Rue et al. [2009] H. Rue, S. Martino, and N. Chopin.  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, 2009.
  * Rue et al. [2022] H. Rue, F. Lindgren, D. Simpson, S. Martino, E. Teixeira Krainski, H. Bakka, A. Riebler, and G.-A. Fuglstad.  _INLA: Full Bayesian Analysis of Latent Gaussian Models using Integrated Nested Laplace Approximations_ , 2022.  R package, version 20.03.17, available from https://www.r-inla.org.
  * Schefzik and Möller [2018] R. Schefzik and A. Möller.  Ensemble postprocessing methods incorporating dependence structures.  In S. Vannitsem, D. S. Wilks, and J. W. Messner, editors, _Statistical Postprocessing of Ensemble Forecasts_ , pages 91–125. Elsevier, 2018.
  * Schefzik et al. [2013] R. Schefzik, T. L. Thorarinsdottir, and T. Gneiting.  Uncertainty quantification in complex simulation models using ensemble copula coupling.  _Statistical science_ , 28(4):616–640, 2013.
  * Schervish [2012] M. J. Schervish.  _Theory of statistics_.  Springer Science & Business Media, 2012.
  * Shores [2007] T. S. Shores.  _Applied linear algebra and matrix analysis_ , volume 2541.  Springer, 2007.
  * South [2017] A. South.  _rnaturalearth: World Map Data from Natural Earth_ , 2017.  URL https://CRAN.R-project.org/package=rnaturalearth.  R package version 0.1.0, Made with Natural Earth naturalearthdata.com.
  * van Straaten et al. [2020] C. van Straaten, K. Whan, D. Coumou, B. van den Hurk, and M. Schmeits.  The influence of aggregation and statistical post-processing on the subseasonal predictability of european temperatures.  _Quarterly Journal of the Royal Meteorological Society_ , 146(731):2654–2670, 2020.
  * Vannitsem et al. [2018] S. Vannitsem, D. S. Wilks, and J. Messner.  _Statistical postprocessing of ensemble forecasts_.  Elsevier, 2018.
  * Vitart et al. [2017] F. Vitart, C. Ardilouze, A. Bonet, A. Brookshaw, M. Chen, C. Codorean, M. Déqué, L. Ferranti, E. Fucile, M. Fuentes, et al.  The subseasonal to seasonal (S2S) prediction project database.  _Bulletin of the American Meteorological Society_ , 98(1):163–173, 2017.
  * Wickham et al. [2019] H. Wickham, M. Averick, J. Bryan, W. Chang, L. D. McGowan, R. François, G. Grolemund, A. Hayes, L. Henry, J. Hester, et al.  Welcome to the Tidyverse.  _Journal of Open Source Software_ , 4(43):1686, 2019.
  * Wilks [2011] D. S. Wilks.  _Statistical methods in the atmospheric sciences_ , volume 100.  Academic press, 2011.
  * Yuen and Stoev [2014] R. Yuen and S. Stoev.  CRPS M-estimation for max-stable models.  _Extremes_ , 17(3):387–410, 2014.

## Supplementary material

##  Appendix A CDS API requests to download reanalysis data

The following python code retrieves the ERA5 temperature reanalysis data via
the Copernicus Data Store (CDS) API444https://cds.climate.copernicus.eu/api-
how-to [Last accessed: 9 August 2022].

    import cdsapi
    c = cdsapi.Client()
    c.retrieve(
      ’reanalysis-era5-single-levels’,
      {
        ’product_type’: ’reanalysis’,
        ’variable’: ’2m_temperature’,
        ’year’: [’2002’, ’2003’, ’2004’, ’2005’, ’2006’,
                 ’2007’, ’2008’, ’2009’, ’2010’, ’2011’,
                 ’2012’, ’2013’, ’2014’, ’2015’, ’2016’,
                 ’2017’, ’2018’, ’2019’, ’2020’, ’2021’],
        ’month’: ’04’,
        ’day’: [’15’, ’16’, ’17’, ’18’, ’19’,
                ’20’, ’21’, ’22’, ’23’, ’24’
        ],
        ’time’: [
            ’00:00’, ’06:00’, ’12:00’, ’18:00’
        ],
        ’format’: ’grib’,
        ’area’: [70, -15, 35, 30],
      },
      ’era5-t2m.grib’)

To retrieve ERA5 total precipitation reanalysis data, change the request as
follows:

    c.retrieve(
      ...
      ’variable’: ’total_precipitation’,
      ...
      ’time’: [’00:00’, ’01:00’, ’02:00’, ’03:00’, ’04:00’, ’05:00’,
               ’06:00’, ’07:00’, ’08:00’, ’09:00’, ’10:00’, ’11:00’,
               ’12:00’, ’13:00’, ’14:00’, ’15:00’, ’16:00’, ’17:00’,
               ’18:00’, ’19:00’, ’20:00’, ’21:00’, ’22:00’, ’23:00’],
      ...
      ’era5-tp.grib’)

##  Appendix B MARS API requests to download ensemble forecasts

The following python code retrieves the 2-metre temperature ensemble forecasts
that was used in this paper from the S2S database via the ECMWF Web
API555https://www.ecmwf.int/en/computing/software/ecmwf-web-api [Last
accessed: 9 August 2022].

    from ecmwfapi import ECMWFDataServer
    server = ECMWFDataServer()
    server.retrieve({
      "class": "s2",
      "dataset": "s2s",
      "date": "2022-04-14",
      "expver": "prod",
      "hdate": "2002-04-14/2003-04-14/2004-04-14/2005-04-14/
                2006-04-14/2007-04-14/2008-04-14/2009-04-14/
                2010-04-14/2011-04-14/2012-04-14/2013-04-14/
                2014-04-14/2015-04-14/2016-04-14/2017-04-14/
                2018-04-14/2019-04-14/2020-04-14/2021-04-14",
      "levtype": "sfc",
      "model": "glob",
      "number": "1/2/3/4/5/6/7/8/9/10",
      "origin": "ecmf",
      "param": "167",
      "step": "0-24/24-48/48-72/72-96/96-120/120-144/
               144-168/168-192/192-216/216-240",
      "stream": "enfh",
      "time": "00:00:00",
      "type": "pf",
      "target": "reforecast-t2m.grib"
    })

To retrieve the precipitation forecasts used in this paper, replace the
following parameters

    server.retrieve({
      ...
      "param": "228228",
      "step": "24/48/72/96/120/144/168/192/216/240",
      "target": "reforecast-tp.grib"
      ...
    })

##  Appendix C Max-and-Smooth R code

We store spatial data in a tibble, ordered by longitude first and latitude
next, with MLEs and observed information matrices stored as nested data in
list columns. The following are the first few lines of the tibble mle used to
store logistic regression parameters and information matrices:

    mle
    # # A tibble: 713  4
    #    long   lat theta_hat J_hat
    #   <dbl> <dbl> <list>    <list>
    # 1   -15  36   <dbl [2]> <dbl [2  2]>
    # 2   -15  37.5 <dbl [2]> <dbl [2  2]>
    # 3   -15  39   <dbl [2]> <dbl [2  2]>
    # #  with 710 more rows

The R code below constructs the sparse precision matrix Q, the sparse
information matrix J_hat, and the vector of MLEs theta_hat. The solve function
provided by the Matrix package is used to solve the sparse linear system to
calculate the smoothed parameter estimates.

    library(Matrix)

    ## joint precision matrix (RW2D block matrix)
    N_i = 23        # number of latitudes (rows)
    N_j = 31        # number of longitudes (columns)
    N   = N_i * N_j # number of grid points
    kappa_hat = c(exp(8.73), exp(6.35)) # optimised hyperparameters
    Di = bandSparse(n=N_i, k=0:1, sym=TRUE,
                    diag=list(c(1, rep(2, N_i-2), 1), rep(-1, N_i-1)))
    Dj = bandSparse(n=N_j, k=0:1, sym=TRUE,
                    diag=list(c(1, rep(2, N_j-2), 1), rep(-1, N_j-1)))
    D  = kronecker(Dj, Diagonal(N_i)) + kronecker(Diagonal(N_j), Di)
    R  = crossprod(D)
    Q  = bdiag(list(kappa_hat[1] * R, kappa_hat[2] * R))

    ## extract vector theta_hat and information block matrix from ‘mle‘
    theta_hat = do.call(c, mle$theta_hat)
    J_hat     = bdiag(mle$J_hat)

    ## theta_hat and J_hat are ordered as (ababab...). we must reorder
    ## them as (a...ab...b) to be compatible with Q
    inds      = c(seq(1, 2*N, 2), seq(2, 2*N, 2))
    theta_hat = theta_hat[inds]
    J_hat     = J_hat[inds, inds]

    ## calculate Max-and-Smooth posterior mean
    theta_ms = drop(solve(Q + J_hat, J_hat %*% theta_hat))

    ## construct parameter vectors that align with the
    ## longitudes and latitudes in ‘mle‘
    alpha_ms = theta_ms[1:N]
    beta_ms  = theta_ms[1:N + N]

##  Appendix D R-INLA smoothing code

We start with a tibble containing latitude/longitude spatial data of MLEs
alpha_hat and their estimation error variances alpha_hat_var.

    alpha_hat
    ## # A tibble: 713  4
    ##    long   lat alpha_hat alpha_hat_var
    ##   <dbl> <dbl>     <dbl>         <dbl>
    ## 1   -15  36       -6.65         15.9
    ## 2   -15  37.5     -3.28          1.43
    ## 3   -15  39       -2.43          1.04
    ## #  with 710 more rows

    ## add the spatial index s used in INLAs RW2D model. the function
    ## as.numeric(factor(...)) assigns increasing integers starting at 1 to latitude
    ## and longitudes.
    N_i = 23
    N_j = 31
    alpha_hat = alpha_hat %>%
      mutate(i = as.numeric(factor(lat)), j = as.numeric(factor(long))) %>%
      mutate(s = (j - 1) * N_i + i)

    ## fit normal measurement error model with RW2D prior. fix error variances using
    ## ‘control.family‘ and ‘scale‘
    inla_alpha = inla(
      formula = alpha_hat ~ -1 + f(s, model=’rw2d’, nrow=N_i, ncol=N_j, constr=FALSE),
      data = alpha_hat,
      family = ’gaussian’,
      control.family = list(hyper=list(prec=list(initial=0, fixed=TRUE))),
      scale = 1 / alpha_hat$alpha_hat_var
    )

    ## extract posterior means of alpha
    alpha_ms = inla_alpha$summary.random$s$mean

The vector alpha_ms contains the smoothed estimates of the intercept
parameters in the same order as the MLEs given in the tibble alpha_hat.

##  Appendix E CRPS and PIT analysis of NGR forecasts

We repeat the leave-one-out verification of temperature forecasts
postprocessed with NGR using the Continuous Ranked Probability Score [CRPS,
Matheson and Winkler, 1976] as the verification metric. For a probabilistic
forecast $f(x)$ issued as a Normal distribution with mean $\mu$ and variance
$\sigma^{2}$, and verifying observation $y$, the CRPS is given by [Gneiting et
al., 2005]

$CRPS(f,y)=\sigma\left\\{z[2\Phi(z)-1]+2\varphi(z)-\pi^{-1/2}\right\\}$ (39)

where $z=(y-\mu)/\sigma$, and $\varphi(x)$ and $\Phi(x)$ are, respectively,
the pdf and cdf of the standard normal distribution.

Figure 10: Spatially averaged CRPS of NGR temperature forecasts. Smaller
values indicate better forecasts.

Figure 10 shows spatially averaged CRPS versus lead time similar to Fig. 8 for
Logscores. In terms of CRPS the independently smoothed postprocessing
parameters with R-INLA perform consistently worse than the the unsmoothed
parameters. When smoothing parameters jointly with Max-and-Smooth, forecasts
improve consistently.

Histograms of Probability Integral Transforms (PITs) are used to assess
forecast calibration [Gneiting et al., 2007]. A forecast at time $t$ and
location $s$ is issued as a cumulative distribution function (cdf)
$F_{s,t}(x)$, and verifying observation $y_{s,t}$ materialised at that
instance. The PIT value of that forecast is then $F_{s,t}(y_{s,t})$. PIT
values of a well-calibrated forecaster have a uniform distribution, and so the
histogram of PIT values should be approximately flat.

Figure 11: PIT histograms of NGR temperature forecasts, separately for each
lead time (rows). Well-calibrated forecasts have flat PIT histograms along the
dashed horizontal line.

Figure 11 shows PIT histograms of (cross-validated) temperature forecasts
postprocessed with NGR. We compare forecasts obtained with unsmoothed
postprocessing parameters (MLE), forecasts obtained with independently
smoothed parameters using R-INLA, and forecasts obtained with parameters
smoothed with Max-and-Smooth. Forecast distributions obtained with unsmoothed
parameters are on average too narrow, as indicated by the overrepresented
outer bins in the PIT histogram. That effect is much reduced in the forecasts
obtained with smoothed parameters. The difference between independently and
dependently smoothed parameters is less strong, but a case can be made that
PIT histograms of Max-and-Smooth forecasts are more uniform than with
indpendent smoothing, especially for lead times greater than 5 days.