# On-chip QNN: Towards Efficient On-Chip Training of
Quantum Neural Networks

1Hanrui Wang∗, 2Zirui Li∗, 3Jiaqi Gu, 4Yongshan Ding, 3David Z. Pan, 1Song Han
1Massachusetts Institute of Technology
2Shanghai Jiao Tong University
3University of Taxes at Austin
4Yale University
https://qmlsys.mit.edu

###### Abstract.

Quantum Neural Network (QNN) is drawing increasing research interest thanks to
its potential to achieve quantum advantage on near-term Noisy Intermediate
Scale Quantum (NISQ) hardware. In order to achieve _scalable_ QNN learning,
the training process needs to be offloaded to real quantum machines instead of
using exponential-cost classical simulators. One common approach to obtain QNN
gradients is _parameter shift_ whose cost scales _linearly_ with the number of
qubits. We present On-chip QNN, the first experimental demonstration of
practical on-chip QNN training with parameter shift. Nevertheless, we find
that due to the significant quantum errors (noises) on real machines,
gradients obtained from naïve parameter shift have low fidelity and thus
degrade the training accuracy. To this end, we further propose _probabilistic
gradient pruning_ to firstly identify gradients with potentially large errors
and then remove them. Specifically, small gradients have larger relative
errors than large ones, thus having a higher probability to be pruned. We
perform extensive experiments on 5 classification tasks with 5 real quantum
machines. The results demonstrate that our on-chip training achieves over 90%
and 60% accuracy for 2-class and 4-class image classification tasks. The
probabilistic gradient pruning brings up to 7% QNN accuracy improvements over
no pruning. Overall, we successfully obtain similar on-chip training accuracy
compared with noise-free simulation but have much better training scalability.
The code for parameter shift on-chip training is available in the TorchQuantum
library.

##  1\. Introduction

∗Equal Contribution.

Quantum Computing (QC) has great potential to achieve exponential acceleration
over classical computers, which represents a computational paradigm shift in
various domains, including chemistry, cryptography, database search, etc.
Among them, Quantum Neural Network (QNN) is one of the most promising
applications to demonstrate quantum advantage on Machine Learning tasks.

In order to achieve QNN quantum advantage, the number of qubit needs to be
large enough, which casts great difficulty in the parameter training process.
In existing QNN work (Harrow et al., 2009; Lloyd et al., 2013), the primary
focus has been building quantum models that can outperform classical model
accuracy. Thus they typically perform training on classical computers through
software simulations and then perform inference with simulators as well
(Figure 1 top). Although classical simulation is useful in understanding the
capabilities of small-size QNN, it is not scalable due to the exponentially
increased time and memory costs ($\mathcal{O}(2^{n}),n$ is the qubit number).
As shown in Figure 2(a), the space (#Regs) and time (#Ops) complexity of
classical simulation grow exponentially as the number of qubits increases. To
the authors’ knowledge, this is the _first experimental demonstration_ of
efficient and scalable QNN on-chip training protocol. The optimization of
parametrized quantum gates is offloaded to the quantum chips with in-situ
gradient computation using _parameter shift_ (Cerezo and others., 2021; Bharti
et al., 2021). We also perform QNN evaluation on real quantum machines, making
the results more _practical_ as in Figure 1 bottom.

Figure 1. In this work, quantum neural network training and inference are both
performed on _real quantum machines_ , making the whole pipeline scalable and
practical.

One of the major challenges to enable scalable and efficient QNN on-chip
learning is the robustness against quantum noise. In the current Noisy
Intermediate Scale Quantum (NISQ) (Preskill, 2018) era, the gate error rates
on real quantum devices are non-negligible ($10^{-3}$ to $10^{-2}$). In the
context of QNN, such errors will lead to _noisy gradients_ which can slow down
convergence or even make training unstable. As shown in Figure 2(b), large
gaps exist between the quantum on-chip training results and the classical
noise-free simulation results.

By carefully investigating the on-chip training process, we observe that small
gradients tend to have large relative variations or even wrong directions
under quantum noises, as shown in Figure 2(c). Also, not all gradient
computations are necessary for the training process, especially for small-
magnitude gradients. Those observations provide great opportunities for us to
boost the robustness and efficiency of QNN on-chip learning. Inspired by that,
we propose a _probabilistic gradient pruning_ method to predict and only
compute gradients of high reliability. Hence we can reduce noise impact and
also save the required number of circuit runs on real quantum machines.

On-chip QNN has following contributions:

  * •

We are the _first work_ to demonstrate the practicality of using parameter
shift on _NISQ machine_ and achieve high accuracy on QNN learning tasks.

  * •

A probabilistic gradient pruning method is put forward to improve the noise
robustness by 5-7% and reduce the #inference on real QC by 2$\times$ while
maintaining the accuracy.

  * •

Experimental deployment on 5 real quantum machines demonstrates that the
proposed method can achieve over 90% and 60% accuracy for 2-class and 4-class
image recognition tasks. Our framework enables scalable, robust, and efficient
training of QNNs with large #qubit and #parameters.

  * •

We integrate the parameter shift on-chip QNN training code to the TorchQuantum
library.

((a))

((b))

((c))

Figure 2. (a) Classical simulation has unscalable computational and memory
costs. (b) Noises create significant accuracy gaps between QNN classical
simulation and on-chip training. (c) Small gradients suffer from larger
relative errors, thus being less reliable.

##  2\. Background

Quantum basics. Quantum circuits use quantum bit (called _qubit_) to store
information, which is a linear combination of two basis states:
$\ket{\psi}=\alpha\ket{0}+\beta\ket{1}$, for $\alpha,\beta\in\mathbb{C}$,
satisfying $|\alpha|^{2}+|\beta|^{2}=1$. An $n$-qubit system can represent a
linear combination of 2n basis states. A 2n-length complex statevector of all
combination coefficients is used to describe the quantum state. To perform
computation on a quantum system, a sequence of parametrized quantum gates are
applied to perform unitary transformation on the statevector, i.e.,
$\ket{\psi(\bm{x},\bm{\theta})}=\cdots
U_{2}(\bm{x},\theta_{2})U_{1}(\bm{x},\theta_{1})\ket{0}$, where $\bm{x}$ is
the input data, and $\bm{\theta}=(\theta_{1},\theta_{2},\dotsc)$ are trainable
parameters in quantum gates. In this way, input data and trainable parameters
are embedded in the quantum state $\ket{\psi(\bm{x},\bm{\theta})}$. The
computation results are obtained by qubit readout which measures the
probability of a qubit state $\ket{\psi}$ collapsing to either $\ket{0}$
(i.e., output $y=+1$) or $\ket{1}$ (i.e., output $y=-1$) according to
$|\alpha|^{2}$ and $|\beta|^{2}$. With sufficient samples, we can compute the
expectation value: $\mathbb{E}[y]=(+1)|\alpha|^{2}+(-1)|\beta|^{2}$. A non-
linear network can be constructed to perform ML tasks by cascading multiple
blocks of quantum gates and measurements.

Quantum noise. In real quantum computer systems, errors (noises) would occur
due to unwanted interactions between qubits, imperfect control signals, or
interference from the environment (Hsieh et al., 2020). For example, quantum
gates introduce _operation errors_ (e.g., coherent errors and stochastic
errors) into the system, and qubits also suffer from _decoherence error_
(spontaneous loss of its stored information) over time. These noisy systems
need to be characterized (Magesan et al., 2012) and calibrated (IBM, [n. d.])
frequently to mitigate the noise impact.

Quantum neural networks. Quantum Machine Learning (QML) (Biamonte et al.,
2017; Wang et al., 2022a, b; Liang et al., 2021; Wang et al., 2021a) aims to
leverage QC techniques to solve machine learning tasks and achieve much higher
efficiency. The path to _quantum advantage_ on QML is typically provided by
the quantum circuit’s ability to generate and estimate highly complex kernels
(Havlíček et al., 2019), which would otherwise be intractable to compute with
conventional computers. They have been shown to have potential speed-up over
classical counterparts in various tasks, including metric learning (Lloyd et
al., 2020), data

Figure 3. Quantum neural network architecture.

analysis (Lloyd et al., 2016). As shown in Figure 3, the quantum neural
network is one type of QML model using variational quantum circuits with
trainable parameters to accomplish feature encoding of input data and perform
complex-valued linear transformations thereafter. Most of QNN trainings are
exploratory and rely on classical simulation of small quantum systems. In our
work, on the contrary, we explore the practical setting: the QNN training and
inference are both performed on real quantum devices.

Pruning. Pruning techniques are widely used in the field of DNN (Han et al.,
2015; Wang et al., 2020b; Wang et al., 2021b; Zhang et al., 2020; Wang et al.,
2020a), performing an important role of the trade-off between accuracy and
memory or time cost (Nguyen-Meidine et al., 2020). Recently, pruning
techniques have been used in quantum tasks. Pruning the ansatz can bring time-
efficient circuit and even higher performance on real QC (Wang et al., 2022a).
In our work, we apply pruning techniques to prune unreliable gradients in
order to mitigate the noise during training.

##  3\. Methodology

To enable QNN on-chip learning, we first introduce an in-situ quantum gradient
computation via parameter shift and its real QC implementation. A
probabilistic gradient pruning method is proposed to save the gradient
computation cost with enhanced noise-robustness and training efficiency.

###  3.1. Parameter shift rule for quantum gradients

Parameter shift rule states that we can calculate the gradient of each
parameter in some quantum circuits by simply shifting the parameter twice and
calculating the difference between two outputs, without changing the structure
of circuits or using any ancilla qubits. Prior works elaborate it based on
quantum circuit function (Crooks, 2019), however, in the next subsection we
will show how parameter shift rules combined with backpropagation can be used
in a real QNN task. Suppose an $m$-qubit quantum circuit is parametrized by
$n$ parameters $\theta=[\theta_{1},\cdots,\theta_{i},\cdots,\theta_{n}]$, the
expectation value of measurements of this circuit can be represented by a
circuit function,

(1) $\small
f(\theta)=\langle\psi|U(\theta_{i})^{{\dagger}}\widehat{Q}U(\theta_{i})|\psi\rangle,\quad
f(\theta)\in\mathbb{R}^{m},\theta\in\mathbb{R}^{n}.$

where $\theta_{i}$ is the scalar parameter whose gradient is to be calculated,
and $U(\theta_{i})$ is the gate where $\theta_{i}$ lies in. Here, for notation
simplicity, we have already absorbed the unitaries before $U(\theta_{i})$ into
$\langle\psi|$, $|\psi\rangle$. Unitaries after $U(\theta_{i})$ and
observables are fused into $\widehat{Q}$. Usually, the gates used in QNN can
be written in the form $U(\theta_{i})=e^{-\frac{i}{2}\theta_{i}H}$. Here $H$
is the Hermitian generator of $U$ with only 2 unique eigenvalues +1 and -1
($H$’s eigenvalues can be $\pm r$, but for simplicity we assume it’s $\pm 1$).
In this way, the gradients of the circuit function $f$ with respect to
$\theta_{i}$ are,

(2) $\displaystyle\frac{\partial
f(\theta)}{\partial\theta_{i}}=\frac{1}{2}\Big{(}f\big{(}\theta_{+}\big{)}-f\big{(}\theta_{-}\big{)}\Big{)},$
$\displaystyle\theta_{+}=[\theta_{1},\cdots,\theta_{i}+\frac{\pi}{2},\cdots,\theta_{n}],\theta_{-}=[\theta_{1},\cdots,\theta_{i}-\frac{\pi}{2},\cdots,\theta_{n}],$

where $\theta_{+}$ and $\theta_{-}$ are the _positive shift_ and _negative
shift_ of $\theta$. Note that this parameter shift rule is _fundamentally
different_ from any numerical difference methods that only approximate the
directional derivatives. Instead, Eq. 2 calculates the _exact_ gradient w.r.t
$\theta_{i}$ without any approximation errors or numerical issues.

We apply softmax on the expectation values of measurements $f(\theta)$ as the
predicted probability for each class. Then we calculate the cross entropy
between the predicted probability distribution $p$ and the target distribution
$t$ as the classification loss $\mathcal{L}$,

(3)
$\small\mathcal{L}(\theta)=-t^{T}\cdot\texttt{softmax}(f(\theta))=-\sum_{j=1}^{m}t_{j}\log{p_{j}},\quad
p_{j}=\frac{e^{f_{j}(\theta)}}{\sum_{j=1}^{m}e^{f_{j}(\theta)}}.$

Then the gradient of the loss function with respect to $\theta_{i}$ is
$\frac{\partial\mathcal{L}(\theta)}{\partial\theta_{i}}=\big{(}\frac{\partial\mathcal{L}(\theta)}{\partial
f(\theta)}\big{)}^{T}\frac{\partial f(\theta)}{\partial\theta_{i}}$.

Here $\frac{\partial f(\theta)}{\partial\theta_{i}}$ can be calculated on real
quantum circuit by the parameter shift rule, and
$\frac{\partial\mathcal{L}(\theta)}{\partial f(\theta)}$ can be efficiently
calculated on classical devices using backpropagation supported by automatic
differentiation frameworks, e.g., PyTorch (Paszke et al., 2019) and TensorFlow
(Abadi et al., 2015).

Now we derive the parameter shift rule used in our QNN models.

Assume $U(\theta_{i})=R_{X}(\theta_{i}),R_{X}(\alpha)=e^{-\frac{i}{2}\alpha
X}$, where $X$ is the Pauli-X matrix.

Firstly, the RX gate is,

(4) $\displaystyle R_{X}(\alpha)$ $\displaystyle=e^{-\frac{i}{2}\alpha
X}=\sum_{k=0}^{\infty}(-i\alpha/2)^{k}X^{k}/k!$
$\displaystyle=\sum_{k=0}^{\infty}(-i\alpha/2)^{2k}X^{2k}/(2k)!+\sum_{k=0}^{\infty}(-i\alpha/2)^{2k+1}X^{2k+1}/(2k+1)!$
$\displaystyle=\sum_{k=0}^{\infty}(-1)^{k}(\alpha/2)^{2k}I/(2k)!-i\sum_{k=0}^{\infty}(-1)^{k}(\alpha/2)^{2k+1}X/(2k+1)!$
$\displaystyle=\cos(\alpha/2)I-i\sin(\alpha/2)X.$

Let $\alpha=\frac{\pi}{2}$, $R_{X}(\pm\frac{\pi}{2})=\frac{1}{\sqrt{2}}(I\mp
iX)$.

As
$f(\theta)=\langle\psi|R_{X}(\theta_{i})^{{\dagger}}\widehat{Q}R_{X}(\theta_{i})|\psi\rangle$,
$R_{X}(\alpha)R_{X}(\beta)=R_{X}(\alpha+\beta)$, and
$\frac{\partial}{\partial\alpha}R_{X}(\alpha)=-\frac{i}{2}XR_{X}(\alpha)$, we
have

(5) $\displaystyle\frac{\partial f(\theta)}{\partial\theta_{i}}=$
$\displaystyle\langle\psi|R_{X}(\theta_{i})^{{\dagger}}(-\frac{i}{2}X)^{{\dagger}}\widehat{Q}R_{X}(\theta_{i})|\psi\rangle+\langle\psi|R_{X}(\theta_{i})^{{\dagger}}\widehat{Q}(-\frac{i}{2}X)R_{X}(\theta_{i})|\psi\rangle$
$\displaystyle=$
$\displaystyle\frac{1}{4}(\langle\psi|R_{X}(\theta_{i})^{{\dagger}}(I-iX)^{{\dagger}}\widehat{Q}(I-iX)R_{X}(\theta_{i})|\psi\rangle$
$\displaystyle-\langle\psi|R_{X}(\theta_{i})^{{\dagger}}(I+iX)^{{\dagger}}\widehat{Q}(I+iX)R_{X}(\theta_{i})|\psi\rangle)$
$\displaystyle=$
$\displaystyle\frac{1}{2}(\langle\psi|R_{X}(\theta_{i})^{{\dagger}}R_{X}(\frac{\pi}{2})^{{\dagger}}\widehat{Q}R_{X}(\frac{\pi}{2})R_{X}(\theta_{i})|\psi\rangle$
$\displaystyle-\langle\psi|R_{X}(\theta_{i})^{{\dagger}}R_{X}(-\frac{\pi}{2})^{{\dagger}}\widehat{Q}R_{X}(-\frac{\pi}{2})R_{X}(\theta_{i})|\psi\rangle)$
$\displaystyle=$ $\displaystyle\frac{1}{2}(f(\theta_{+})-f(\theta_{-})).$

Without loss of generality, the derivation holds for all unitaries of the form
$e^{-\frac{i}{2}\alpha H}$, e.g., RX, RY, RZ, XX, YY, ZZ, where $H$ is a
Hermitian matrix with only 2 unique eigenvalues +1 and -1.

In our circuit functions, we assume each parameter lies in exactly one gate.
However, there are cases that one parameter lies in multiple gates. In that
case, we only need to calculate the gradient of the parameter in those gates
separately and sum the gradients up to get the gradient of that parameter.

Figure 4. Quantum gradient calculation using parameter shift on real quantum
devices.

###  3.2. In-situ gradient computation on real QC

To realize QNN on-chip learning, we implement a TrainingEngine, described in
Alg. 1. This TrainingEngine contains three parts.

Jacobian calculation via parameter shift. In the first part, we sample a mini-
batch of training data $\mathcal{I}$ in Line 6. For each example of the mini-
batch, we set up the quantum encoder gates and then iteratively evaluate
gradients for all parameters. In each iteration, we shift the parameter
$\theta_{i}$ twice by $+\pi/2$ and $-\pi/2$ respectively. After each shift, we
execute the shifted circuit on quantum hardware. The circuit will be _created,
validated, queued, and finally run on real quantum machines_. As soon as we
get the returned results of the two shifted circuits, i.e., $f(\theta_{+})$
and $f(\theta_{-})$, we apply Eq. 2 to obtain the upstream gradient
$\frac{\partial f(\theta)}{\partial\theta_{i}}$, illustrated in the left part
of Figure 4. Finally, we obtain the Jacobian matrix $\frac{\partial
f(\theta)}{\partial\theta}$.

Figure 5. Efficient on-chip quantum gradient calculation with probabilistic
gradient pruning. Gradient magnitudes are accumulated within the accumulation
window and used as the sampling distribution. Based on the distribution,
gradients are probabilistically pruned with a ratio $r$ in the pruning window
to mitigate noises and stabilize training.

Down-stream gradient backpropagation. In the second part, we run the circuit
without shift and get the measurement result $f(\theta)$. Then we apply
softmax and cross-entropy function to the measured logits. In the end, we get
the training loss $\mathcal{L}(\theta)$. Then we run back-propagation _only
from the loss to the logits_ to get the down-stream gradients
$\frac{\partial\mathcal{L}(\theta)}{\partial f(\theta)}$, shown in the right
part of Figure 4.

Gradient calculation. In the third part, we calculate the dot-product between
down-stream gradients and the Jacobian and get the final gradients
$\frac{\partial\mathcal{L}(\theta)}{\partial\theta}=(\frac{\partial
f(\theta)}{\partial\theta})^{T}\frac{\partial\mathcal{L}(\theta)}{\partial
f(\theta)}$.

###  3.3. Probabilistic quantum gradient pruning

On quantum chips, there exist various noises and errors that could potentially
diminish the fidelity of the computation results. When the gradient magnitude
is small, noises could easily overwhelm the signals, such that the _gradients
calculated on real quantum circuit become unreliable when they have small
magnitude_. Those unreliable gradients have harmful effects on training
convergence. Skipping the evaluation on those unreliable gradients can benefit
both training convergence and efficiency. Besides, we observe that for most
parameters, if the gradient magnitudes are far from zero for several steps, it
will likely keep far from zero in the next several steps. Similarly, if the
gradient magnitude remains small for some steps, it will likely keep small in
the next several steps. This means the _gradient reliability is predictable_
to some extent. Therefore, we propose the gradient pruning method to sample
the parameters whose gradients are more reliable. This method helps training
converge faster while also saving time by skipping the evaluation of
unreliable gradients.

Alg. 1 describes the QNN on-chip training flow with probabilistic gradient
pruning. We divide all the training steps into $S$ stages and perform the
pruning method periodically on each stage. For every stage, we split it into
two phases, shown in Figure 5. The first phase is called _magnitude
accumulation_ with an accumulation window width $w_{a}$, and the second is
called _probabilistic gradient pruning_ (PGP) with a pruning window width
$w_{p}$. We only apply pruning in the second phase, while the parameter subset
is sampled from a probability distribution $\tilde{\theta}=\\{\theta_{i}\sim
P_{M}(\theta)|1\leq i\leq(1-r)n\\}$ based on the gradient information
collected within the accumulation window.

In Lines 4-9, within the accumulation window, we record the magnitude of
gradients of each parameter in each step and accumulate them until the window
is over. At the end of the first phase, we can get an accumulator $M$ that
records the accumulated gradient magnitude for each parameter. Thus, when the
pruning phase starts, we normalize the accumulated gradient magnitude and pass
it to our sampler as the sampling distribution. In each pruning step, the
sampler samples a subset of parameters $\tilde{\theta}$ with a pruning ratio
of $r$, and we only evaluate gradients for them while the rest
$\theta\backslash\tilde{\theta}$ is temporarily frozen.

There are three important hyper-parameters in our gradient pruning method: 1)
accumulation window width $w_{a}$, 2) pruning ratio $r$, and 3) pruning window
width $w_{p}$. The accumulation window width and pruning window width decide
the reliability of the gradient trend evaluation and our confidence in it,
respectively. The pruning ratio can be tuned to balance the gradient variances
caused by noise perturbation and pruning. Thus, the percentage of the time
saved by our probabilistic gradient pruning method is
$r\frac{w_{p}}{w_{a}+w_{p}}\times 100\%$. In our experiments, we find that the
setting ($w_{a}$=1, $w_{p}$=2$\sim$3, $r$=0.3$\sim$0.5) usually works well in
all cases.

Input : Accumulation window width $w_{a}$, gradient pruning ratio $r$, pruning
window width $w_{p}$, training objective $\mathcal{L}$, initial parameters
$\theta^{0}\in\mathbb{R}^{n}$, training data $\mathcal{D}_{trn}$, initial step
size $\eta^{0}$, and total stages $S$.

$\theta\leftarrow\theta^{0},~{}~{}\eta\leftarrow\eta^{0}~{}~{}t\leftarrow 0$;

for  _$s=1,2,.\cdots,S$_ do

Initialize gradient magnitude accumulator $M\leftarrow\mathbf{0}^{n}$;

for  _$\tau_{a}=1,2,\cdots,w_{a}$_ do

$t\leftarrow t+1$;

Sample a mini-batch $\mathcal{I}\sim\mathcal{D}_{trn}$;

In-situ gradient evaluation via parameter shift
$\nabla_{\theta}\mathcal{L}_{\mathcal{I}}(\theta)=\frac{1}{2}(\frac{\partial
f(\theta)}{\partial\theta})^{T}\frac{\partial\mathcal{L}(\theta)}{f(\theta)}$;

Parameter update:
$\theta\leftarrow\theta-\eta\nabla_{\theta}\mathcal{L}_{\mathcal{I}}(\theta)$;

Update magnitude accumulator $M\leftarrow
M+|\nabla_{\theta}\mathcal{L}_{\mathcal{I}}(\theta)|$;

for  _$\tau_{p}\leftarrow 1,2,\cdots,w_{p}$_ do

$t\leftarrow t+1$;

Sample a mini-batch $\mathcal{I}\sim\mathcal{D}_{trn}$;

Sample a subset with a ratio $r$ based on accumulated gradient magnitude:
$\tilde{\theta}=\\{\theta_{i}\sim P_{M}(\theta)|1\leq i\leq(1-r)n\\}$;

$\tilde{\theta}\leftarrow\tilde{\theta}-\eta\nabla_{\tilde{\theta}}\mathcal{L}_{\mathcal{I}}(\theta)$;

Output : Converged parameters $\theta$

Algorithm 1 QNN On-Chip Training with Probabilistic Gradient Pruning

##  4\. Experiments

In this section, we deploy our QNN on-chip learning framework on real QC and
evaluate it on 5 QML tasks for image and vowel recognition. Compared with
classical QNN training protocols, we can achieve 2-4% real QC test accuracy
improvement with 2$\times$ convergence speedup. We also conduct extensive
ablation studies to validate our scalability and the effectiveness of the
proposed probabilistic gradient pruning method.

###  4.1. Experiment Setups

Benchmarks. We conduct our experiments on 5 QML tasks. QML are all
classification tasks including MNIST (Lecun et al., 1998) 4-class (0, 1, 2,
3), 2-class (3 and 6); Fashion (Xiao et al., 2017) 4-class (t-shirt/top,
trouser, pullover, dress), 2-class (dress and shirt); Vowel (Deterding, 1989)
4-class(hid, hId, had, hOd). MNIST and Fashion 2-class use the front 500
images as the training set and randomly sampled 300 images as the validation
set. MNIST, Fashion 4-class uses the front 100 images as the training set and
also randomly sampled 300 images as the validation set. The input images are
all $28\times 28$. We firstly center-crop them to $24\times 24$ and then down-
sample them to $4\times 4$ for MNIST and Fashion 2 and 4-class tasks. Vowel
4-class uses the front 100 samples as the training set and randomly sampled
300 samples as the validation set. For each sample, we perform principal
component analysis (PCA) for the vowel features and take the 10 most
significant dimensions.

Table 1. Accuracy comparison among different settings. Sim. represents
simulation.

Method | Acc. | MNIST-4 | MNIST-2 | Fashion-4 | Fashion-2 | Vowel-4
---|---|---|---|---|---|---
|  | Jarkata | Jarkata | Manila | Santiago | Lima
Classical-Train | Simu. | 0.61 | 0.88 | 0.73 | 0.89 | 0.37
Classical-Train |  | 0.59 | 0.79 | 0.54 | 0.89 | 0.31
QC-Train | QC | 0.59 | 0.83 | 0.49 | 0.84 | 0.34
QC-Train-PGP |  | 0.64 | 0.86 | 0.57 | 0.91 | 0.36

((a))

((b))

Figure 6. Real QC validation accuracy curves on different datasets and
different quantum devices.

All the tasks use four logical qubits. To embed classical image and vowel
features to the quantum states, we first flatten them and then encode them
with rotation gates. For down-sampled $4\times 4$ images, we use 4RY, 4RZ,
4RX, and 4RY gates as the encoder. We put the 16 classical input values to the
phases of 16 rotation gates, respectively. Therefore we can encode the
classical values to quantum states. For 10 vowel features, we use 4RY, 4RZ,
and 2RX gates for encoding.

The encoding gates are our hand-designed circuits. Our circuits are composed
of several layers. There are 7 kinds of layers used to construct our circuits.
($i$) RX layer: Add RX gates to all wires; ($ii$) RY layer: same structure as
in RX layer; ($iii$) RZ layer: same structure as in RX layer; ($iv$) RZZ
layer: add RZZ gates to all logical adjacent wires and the logical farthest
wires to form a ring connection, for example, an RZZ layer in a 4-qubit
circuit contains 4 RZZ gates which lie on wires 1 and 2, 2 and 3, 3 and 4, 4
and 1; ($v$) RXX layer: same structure as in RZZ layer; ($vi$) RZX layer: same
structure as in RZZ layer; ($vii$) CZ layer: add CZ gates to all logical
adjacent wires.

For MNIST and Fashion 2-class tasks, the circuit contains 1 RZZ layer followed
by 1 RY layer. For MNIST 4-class task, the circuit contains 3 RX+RY+RZ+CZ
layers (1 RX layer, 1 RY layer, 1 RZ layer, and 1 CZ layer in series). For
Fashion 4-class task, the circuit contains 3 RZZ+RY layers (1 RZZ layer
followed by 1 RY layer). For Vowel 4-class task, the circuit contains 2
RZZ+RXX layers (1 RZZ layer followed by 1 RXX layer).

For the output of our quantum circuits, we measure the expectation values on
Pauli-Z basis and obtain a value [-1, 1] from each qubit. For 2-class, we sum
the qubit 0 and 1, 2, and 3 respectively to get 2 output values. For 4-class,
we just use the four expectation values as 4 output values. Then we process
the output values by Softmax to get probabilities.

Quantum devices and compiler configurations. We use IBM quantum computers via
qiskit API (IBM, [n. d.]) to submit our circuits to real superconducting
quantum devices and achieve quantum on-chip training. We set all the circuits
to run 1024 shots.

Baseline. We have two baselines. (1) QC-Train: We train our model without
gradient pruning, i.e., calculating gradients of every parameter in each step.
The gradient calculation is deployed on real quantum circuits. (2) Classical-
Train: We train our QNN model completely on classical computers. We use a
vector to record the amplitudes of the quantum state, utilize complex matrix
multiplication to simulate quantum gates, and sample based on the amplitude
vector to simulate quantum measurement.

The QC-Train-PGP line shows training on real quantum circuits while applying
our probabilistic gradient pruning. In all the cases, we adopt accumulation
window size 1, pruning ratio 0.5, and pruning window size 2, except for
Fashion-4, we adopt pruning ratio 0.7, and other settings remain the same.

###  4.2. Main Results

QNN results. Table 1 shows the accuracy of comparison on 5 tasks. In each
task, we show 4 accuracy values, which are (1) accuracy of Classical-Train
tested on classical devices, (2) accuracy of Classical-Train tested on real
quantum circuits; (3) accuracy of QC-Train tested on real quantum circuits;
(4) accuracy of QC-Train-PGP tested on real quantum circuits. In each task,
the accuracy is collected after finishing a certain number of circuit runs. We
train and evaluate MNIST-2 and MNIST-2 on ibmq_jakarta, Fashion-4 on
ibmq_manila, Fashion-2 on ibmq_santiago, and Vowel-4 on ibmq_lima.

The noise-free accuracy is usually the highest among the other 3 accuracy,
because it represents the accuracy without any noise perturbation. The QC-
Train-PGP usually takes second place because compared to Classical-Train, it
has the advantage of noise awareness, and compared to QC-Train, it suffers
less from noise thanks to gradient pruning.

Figure 7. Ablation on pruning ratio, accumulation window width, and pruning
window width.

Training curves. Figure 6 shows the real QC validation accuracy curve during
training. The X-axis is the number of inferences (how many circuits have been
run). The Y-axis is the accuracy of the validation dataset tested on real
quantum circuits. MNIST 4-class runs on the ibmq_jakarta machine. We observe
that given a fixed inference budget, our QC-Train-PGP achieves the best
accuracy of 63.7% while the Classical-Train only achieves 59.3%.

We further train Fashion 2-class on ibmq_santiago. QC-Train-PGP only takes
13.9k inferences to reach the peak accuracy 90.7%, while the best accuracy
Classical-Train can achieve is merely 88.7% at the cost of over 30k
inferences.

###  4.3. Ablation Studies

Ablation on gradient pruning. In Figure 7, we evaluate the training
performance with different pruning ratios $r$, accumulation window size
$w_{a}$, and pruning window size $w_{p}$ on Fashion-4 and MNIST-2 tasks. We
find that the $r=0.5$ is generally a good setting for our tasks. Overly large
pruning ratios will induce too many gradient variances that harm the training
convergence. For the accumulation window size, $w_{a}=1$ or $2$ are suitable
choices. When $w_{a}$ is too large, the accumulated gradient magnitudes are
similar among all parameters, leading to a nearly uniform sampling
distribution. This will bring undifferentiated pruning, and the accuracy will
drop as the Fashion-4 curve shows. The pruning window $w_{p}$ should also not
be too large. As $w_{p}$ grows, the accumulated gradient magnitudes used to
instruct our pruning become less reliable.

Discussion on scalability. Figure 8 shows the superior scalability of quantum
on-chip training. Classical simulation runtime exponentially increases as
#qubits scales up,

Table 2. Proposed probabilistic pruning is better than deterministic pruning.

Method | MNIST-4 | MNIST-2 | Fashion-4 | Fashion-2
---|---|---|---|---
Deterministic | 0.61 | 0.82 | 0.72 | 0.89
Probabilistic | 0.62 | 0.85 | 0.79 | 0.90

Table 3. Adam optimizer can outperform SGD and Momentum optimizers.

Optimizer | MNIST-4 | MNIST-2 | Fashion-4 | Fashion-2
---|---|---|---|---
SGD | 0.5 | 0.8 | 0.45 | 76
Momentum | 0.55 | 0.83 | 0.66 | 0.90
Adam | 0.61 | 0.88 | 0.75 | 0.91

while the runtime on real quantum machines scales nearly linearly to #qubits.
The classical curve in Figure 8 represents runtime and memory cost of running
50 circuits of different #qubits with 16 rotation gates and 32 RZZ gates. The
curve before 22 qubits is measured on a single NVIDIA RTX 2080 Ti GPU; points
after 24 qubits are extrapolated. The quantum curve before 27 qubits is tested
on ibmq_toronto; the points after 30 qubits are extrapolated.

We can observe clear quantum advantages on circuits with more than 27 qubits.
In terms of memory cost, classical simulation consumes thousands of Gigabits
for storage which is intractable. In contrast, on quantum machines, the
information is stored in the quantum state of the circuit itself with
negligible memory cost.

Probabilistic vs. deterministic gradient pruning. Our pruning is decided by a
random sampler based on the accumulated gradient magnitude. We call this
probabilistic pruning. If the sampler only samples the parameters with the
biggest accumulated gradient magnitude, this is called deterministic pruning.
We adopt probabilistic pruning instead of deterministic pruning because
deterministic pruning limits the degree of freedom and increases the gradient
sampling bias. Table 4.3 shows that deterministic pruning has 1%-7% accuracy
loss compared with probabilistic pruning.

Different optimizers. Table 3 shows the accuracy tested on classical devices
trained with different optimizers. The learning rate is controlled by a cosine
scheduler from 0.3 in the beginning to 0.03 in the end. We test SGD, SGD with
a momentum factor of 0.8 (Qian, 1999), and Adam on MNIST-4, MNIST-2,
Fashion-4, and Fashion-2, and found that Adam always performs the best. Hence,
all the experiments are done using Adam (Kingma and Ba, 2015) optimizers by
default.

Figure 8. Runtime (s) and memory cost comparison between classical simulation
and quantum on-chip run.

##  5\. Conclusion

In this work, for the first time, we present an efficient and robust on-chip
training framework for quantum neural networks (QNN) and demonstrate its
effectiveness on real quantum devices. By leveraging parameter shift, we can
calculate the exact quantum gradients directly on quantum machines, thus
achieving high _scalability_. To alleviate the negative impact of quantum
noises on gradients, we further propose the probabilistic gradient pruning
technique to avoid updating parameters with unreliable gradients. Experimental
results on 5 classification tasks and 5 machines demonstrate that our method
can achieve comparable accuracy with pure noise-free simulation. We hope this
work can open an avenue towards practical training of large QNN models for
quantum advantage.

## References

  * (1)
  * Abadi et al. (2015) Martín Abadi et al. 2015\.  TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems.  Software available from tensorflow.org.
  * Bharti et al. (2021) Kishor Bharti et al. 2021\.  Noisy intermediate-scale quantum (NISQ) algorithms.  arXiv:2101.08448 [quant-ph]
  * Biamonte et al. (2017) Jacob Biamonte, Peter Wittek, Nicola Pancotti, Patrick Rebentrost, Nathan Wiebe, and Seth Lloyd. 2017.  Quantum machine learning.  _Nature_ 549, 7671 (2017).
  * Cerezo and others. (2021) M. Cerezo and others. 2021\.  Variational quantum algorithms.  _Nature Reviews Physics_ 3, 9 (Aug 2021), 625–644.  https://doi.org/10.1038/s42254-021-00348-9
  * Crooks (2019) Gavin E Crooks. 2019\.  Gradients of parameterized quantum gates using the parameter-shift rule and gate decomposition.  _arXiv:1905.13311_ (2019).
  * Deterding (1989) D.H. Deterding. 1989\.  Speaker Normalisation for Automatic Speech Recognition.  PhD thesis, University of Cambridge.
  * Han et al. (2015) Song Han, Huizi Mao, and William J Dally. 2015.  Deep compression: Compressing deep neural networks with pruning, trained quantization and huffman coding.  _arXiv preprint arXiv:1510.00149_ (2015).
  * Harrow et al. (2009) Aram W Harrow, Avinatan Hassidim, and Seth Lloyd. 2009\.  Quantum algorithm for linear systems of equations.  _Physical review letters_ 103, 15 (2009), 150502\.
  * Havlíček et al. (2019) Vojtěch Havlíček et al. 2019\.  Supervised learning with quantum-enhanced feature spaces.  _Nature_ 567, 7747 (2019), 209–212.
  * Hsieh et al. (2020) Cheng-Yun Hsieh, Chen-Hung Wu, Chia-Hsien Huang, His-Sheng Goan, and James Chien Mo Li. 2020.  Realistic fault models and fault simulation for quantum dot quantum circuits. In _2020 57th (DAC)_. IEEE, 1–6.
  * IBM ([n. d.]) Qiskit IBM. [n. d.].
  * Kingma and Ba (2015) D. Kingma and J. Ba. 2015.  Adam: A Method for Stochastic Optimization. In _Proceedings of the International Conference on Learning Representations_.
  * Lecun et al. (1998) Y. Lecun, L. Bottou, Y. Bengio, and P. Haffner. 1998\.  Gradient-based learning applied to document recognition.  _Proc. IEEE_ 86, 11 (1998), 2278–2324.
  * Liang et al. (2021) Zhiding Liang, Zhepeng Wang, Junhuan Yang, Lei Yang, Yiyu Shi, and Weiwen Jiang. 2021\.  Can Noise on Qubits Be Learned in Quantum Neural Network? A Case Study on QuantumFlow. In _2021 IEEE/ACM International Conference On Computer Aided Design (ICCAD)_. IEEE, 1–7.
  * Lloyd et al. (2016) Seth Lloyd, Silvano Garnerone, and Paolo Zanardi. 2016\.  Quantum algorithms for topological and geometric analysis of data.  _Nature communications_ 7, 1 (2016).
  * Lloyd et al. (2013) Seth Lloyd, Masoud Mohseni, and Patrick Rebentrost. 2013\.  Quantum algorithms for supervised and unsupervised machine learning.  _arXiv:1307.0411_ (2013).
  * Lloyd et al. (2020) Seth Lloyd, Maria Schuld, Aroosa Ijaz, Josh Izaac, and Nathan Killoran. 2020.  Quantum embeddings for machine learning.  _arXiv:2001.03622_ (2020).
  * Magesan et al. (2012) Easwar Magesan, Jay M Gambetta, and Joseph Emerson. 2012\.  Characterizing quantum gates via randomized benchmarking.  _Physical Review A_ 85, 4 (2012).
  * Nguyen-Meidine et al. (2020) Le Thanh Nguyen-Meidine et al. 2020\.  Progressive Gradient Pruning for Classification, Detection and DomainAdaptation.  arXiv:1906.08746 [cs.LG]
  * Paszke et al. (2019) Adam Paszke et al. 2019\.  Pytorch: An imperative style, high-performance deep learning library.  _arXiv preprint arXiv:1912.01703_ (2019).
  * Preskill (2018) John Preskill. 2018\.  Quantum Computing in the NISQ era and beyond.  _Quantum_ 2 (2018), 79.
  * Qian (1999) Ning Qian. 1999\.  On the momentum term in gradient descent learning algorithms.  _Neural networks_ (1999).
  * Wang et al. (2020a) Hanrui Wang et al. 2020a.  Efficient algorithms and hardware for natural language processing.  (2020).
  * Wang et al. (2022a) Hanrui Wang, Yongshan Ding, Jiaqi Gu, Yujun Lin, David Z Pan, Frederic T Chong, and Song Han. 2022a.  Quantumnas: Noise-adaptive search for robust quantum circuits.  _2022 IEEE International Symposium on High-Performance Computer Architecture (HPCA)_ (2022).
  * Wang et al. (2022b) Hanrui Wang, Jiaqi Gu, Yongshan Ding, Zirui Li, Frederic T Chong, David Z Pan, and Song Han. 2022b.  RobustQuantum: Noise-Aware Training for Robust Quantum Neural Networks.  _The 59th Design Automation Conference (DAC 2022)_ (2022).
  * Wang et al. (2021b) Hanrui Wang, Zhekai Zhang, and Song Han. 2021b.  Spatten: Efficient sparse attention architecture with cascade token and head pruning. In _2021 IEEE International Symposium on High-Performance Computer Architecture (HPCA)_. IEEE, 97–110.
  * Wang et al. (2020b) Tianzhe Wang, Kuan Wang, Han Cai, Ji Lin, Zhijian Liu, Hanrui Wang, Yujun Lin, and Song Han. 2020b.  Apq: Joint search for network architecture, pruning and quantization policy. In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_. 2078–2087.
  * Wang et al. (2021a) Zhepeng Wang, Zhiding Liang, Shanglin Zhou, Caiwen Ding, Yiyu Shi, and Weiwen Jiang. 2021a.  Exploration of Quantum Neural Architecture by Mixing Quantum Neuron Designs. In _2021 IEEE/ACM International Conference On Computer Aided Design (ICCAD)_. IEEE, 1–7.
  * Xiao et al. (2017) Han Xiao, Kashif Rasul, and Roland Vollgraf. 2017.  Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms.  _arXiv:1708.07747_ (2017).
  * Zhang et al. (2020) Zhekai Zhang, Hanrui Wang, Song Han, and William J Dally. 2020\.  Sparch: Efficient architecture for sparse matrix multiplication. In _2020 IEEE International Symposium on High Performance Computer Architecture (HPCA)_. IEEE, 261–274.


# Horseshoe Priors for Time-Varying AR and GARCH Processes

John W. G. Addy  Intelligent Data Ecosystems, Rothamsted Research, UK  Chloe
Maclaren  Department of Crop Production Ecology, Swedish University of
Agricultural Science, Sweden  Kirsty Hassall  Intelligent Data Ecosystems,
Rothamsted Research, UK

###### Abstract

Grassland ecosystems support a wide range of species and provide key services
including food production, carbon storage, biodiversity support, and flood
mitigation. However, yield stability in these grassland systems is not yet
well understood, with recent evidence suggesting water stress throughout
summer and warmer temperatures in late summer reduce yield stability. In this
study we investigate how grassland yield stability of the Park Grass
Experiment, UK, has changed over time by developing a Bayesian time-varying
Autoregressive and time-varying Generalised Autoregressive Conditional
Heterogeneity model using the variance-parameterised Gamma likelihood
function.

##  1 Introduction

Climate change is pushing Earth’s biosphere into conditions that have not been
experienced for millennia and that are novel to present-day ecosystems (Bova
et al., 2021). To forecast how ecosystems will respond, a robust understanding
of their stability is needed. Less stable ecosystems are more prone to
disruptions in or losses of function, posing a greater risk to the
biodiversity and human communities that depend on those ecosystems for life-
supporting goods and services (Willcock et al., 2023). Grasslands are a
globally important ecosystem that support a wide range of species and provide
key services to human society including food production, carbon storage,
biodiversity support, and flood mitigation (Bengtsson et al., 2019). In the
UK, grasslands cover nearly 40% of the land area, and may be highly vulnerable
to climate change. Grassland productivity in the south-east of England
forecast to decline by as much as 48-50% by 2080 if global temperature
increases exceed 1.5∘C (Addy et al., 2022).

Previous studies have found that lower grassland yields in the south-east of
England are associated with warmer and drier years, which are predicted to
increase under climate change (Addy et al., 2021; Cashen, 1947). However, the
yield stability (or variance in yield over time) in these grassland systems is
not yet well understood. Recent empirical evidence suggests that water stress
throughout summer (between March and October) and warmer temperatures in late
summer reduce yield stability (Macholdt et al., 2023), but these trends have
not yet been formally modelled through time-series processes.

Yields from grasslands ecosystems are known to follow an autoregressive lag-1
process and are highly correlated with yield in the previous year (Jenkinson
et al., 1994). This is likely due to high-yielding years resulting in higher
seed production or higher allocation of resources to energy-storage organs
(such as rhizomes), allowing the grassland to make better use of the
conditions and available resources in the following year (Ma et al., 2021).
Such processes may also lead to autoregression in yield stability. Seeds
produced or energy stored in good years could provide a buffer to unfavourable
conditions in subsequent years, while low seed production and scarce energy
storage in bad years could mean grasslands are subsequently more susceptible
to the effects of favourable vs poor conditions. Similarly, shifts in
grassland community composition and diversity in response to changing
conditions may affect the overall productivity and/or resilience of the
ecosystem, leading to patterns in stability over time (Zhang et al., 2022,
2018).

In this study, we seek greater insight into how grassland yield stability has
responded to change over time by using Autoregressive Conditional
Heterogeneity (ARCH) and Generalised Autoregressive Conditional Heterogeneity
(GARCH) processes. These are further extensions of time-series models enabling
a formal approach to modelling yield stability. An ARCH process occurs when
the squared residual error $\varepsilon^{2}$ is autoregressive (Engle, 1982).
In addition to an ARCH process, a GARCH process also considers an additional
moving average process of the error variance $\sigma^{2}$ (Bollerslev, 1986).
Time-Varying Autoregression (TV-AR) and Time-Varying Generalised
Autoregressive Conditional Heterogeneity (TV-GARCH) models estimate time-
dependent coefficients similar to time-varying regression models (Zhang & Wu,
2012). Although there are methods which consider the yield stability of
cropping systems (Piepho, 1998), these methods do not consider the time-
varying GARCH process of grassland yields. Incorporating a TV-AR and a TV-
GARCH process can provide insights into how sensitive mean and variance
estimates, and therefore yield stability of grassland yields, are over time.

With TV-AR and TV-GARCH processes being powerful functions which assess how AR
and GARCH terms change over time, they are typically estimated by smoothed
non-parametric functions (Hoover et al., 1998) which may need regulating to
prevent over-fitting of the data (Tibshirani, 1996). If smoothed non-
parametric functions are not regulated this can lead to unrealistic estimates
and the misinterpretation of the smoothed relationship. These smoothed non-
parametric functions can be considered as a linear combination of basis
functions and basis coefficients, penalised-splines (P-splines) have
introduced the global shrinkage of basis coefficients as a method to avoid
over-fitting (Eilers & Marx, 1996). Bayesian approaches consider priors for
P-splines to assume a zero-centred multivariate-Normal distribution on spline
basis coefficients (Wood, 2017) with a user defined penalty matrix and global
smoothing parameter representing the covariance structure (Lang & Brezger,
2004). However, the method of P-splines and global shrinkage can dominate the
smoothing of a non-parametric functions. Horseshoe priors (Carvalho et al.,
2010) have been developed to allow more localised shrinkage for each
coefficient through a series of zero-centred Normal distributions. Although
they are more robust in dealing with over-fitting, posterior estimates from
Horseshoe priors are often correlated with no formalised prior structure and
little understanding of the influence of covariate shrinkage introduced by
localised shrinkage parameters. We compare how different prior specifications
regulate smoothed non-parametric functions for TV-AR and TV-GARCH processes.

To compare how different prior specifications regulate smoothed non-parametric
functions, we first generalise the global shrinkage parameter as a covariance
matrix by introducing an Inverse-Wishart prior (Kass & Natarajan, 2006) on the
variance component of a zero-centred multivariate-Normal distribution.
Additionally, we develop and apply multivariate Horseshoe priors for smooth
non-parametric time-varying functions by introducing the LKJ-prior
(Lewandowski et al., 2009) on the correlation structure of a decomposed
covariance matrix (Barnard et al., 2000). The influence on the effect of
global and localised shrinkage is compared for these distinct prior
specifications. Two simulation studies are presented addressing how different
prior specifications for basis coefficients of the non-parametric TV-AR(1) and
TV-GARCH(1, 1) processes influence model prediction and regulate the estimated
time-varying function.

With TV-AR and TV-GARCH processes used in economics and the medical sciences
(Giudici et al., 2023), we use these processes to better understand and
formalise how autoregressive and conditional heterogeneous terms change over
time for grassland yield using data from the Park Grass Experiment,
Hertfordshire. To ensure mean estimates of yield are strictly non-negative, we
develop a Bayesian TV-AR and TV-GARCH model using the variance-parameterised
Gamma likelihood function (Addy et al., 2022) along with a mean-variance
relationship across various treatments of the Park Grass Experiment.

##  2 Data

###  2.1 Park Grass Experiment

The herbage yield data analysed in this work comes from the Park Grass
Experiment at Rothamsted Research. Different fertilizer inputs affect the
herbage yield from each plot (Figure 1) along with the species composition of
the original meadowland, with different species composition across plots
(Macdonald, 2018). Herbage yield data were taken from five treatments, two
with no fertilizer inputs (plots 2.2 and 3), nitrogen inputs at 96 kg N ha -1
plus minerals (plot 14.2), 48 kg N ha -1 plus minerals (plot 16) and mineral
only plot (plot 7). Each treatment has four liming sections with a maintained
pH of 7, 6, 5 and unlimed. There are 20 sections of Park Grass used in this
analysis (five treatments and 4 liming sections) with data from 1966 to 2018.
Data from 2003 was missing for plots 14.2, 16, 2.2, 7 that were maintained at
pH levels 6 and unlimed. Also from 2003 data was missing for plot 3 unlimed.
Missing data was imputed using mean estimates from regression terms and TV-
AR(1) estimates (Little & Rubin, 1986).

One of the main objectives of this study is to apply time-varying time-series
approaches to herbage yield data to better understand, integrate and formally
model how AR and GARCH processes have changed over time and if these processes
have been influenced by climate or environmental change. For example, yields
from the late-1990s to early-2000s were slightly more variable compared to any
other time-period (Figure 1). We also develop a TV-GARCH model including a
heteroskedastic process fitted as a mean-variance relationship with a Gamma
model. It is known that climatic conditions affect grassland production and
with global temperatures likely to exceed 1.5∘C by 2100 (IPCC, 2023),
understanding how AR and GARCH processes change with time will better inform
how grassland production will change over time. Seasonal mean temperature (∘C)
and total rainfall (mm) were provided by the Rothamsted Meteorological Station
from 1966 to 2018.

##  3 Smoothed Functions and Priors

###  3.1 Non-Parametric Smoothed Functions

Consider $f(x)$ to be a smooth non-parametric function between predictor $x$
and response $y$. Estimates of the smoothed function $f(x)$ can be considered
as a linear combination of basis functions $B_{p}$ and basis coefficients
$\theta_{p}$, over $p$ basis dimensions,

$f(x_{i}|\theta)=\sum_{p=1}^{m}\theta_{p}B_{p}(x_{i})$

(Wood, 2017). In this study we only consider $p$ equidistant Gaussian
functions over a total of $m$ dimensions. Although the method of non-
parametric smoothing is considered flexible, when $\theta_{p}$ is un-penalised
this can result in the smoothed function $f(x)$ over-fitting the data.
Consequences of over-fitting can result in unrealistic estimates of the smooth
relationship which influences the interpretation of the model and the overall
conclusions of the analysis.

###  3.2 P-Splines Priors and Global Shrinkage

Shrinking of the basis coefficient’s sample space compensates for model over-
fitting by narrowing the plausible sample region for each coefficient. This
results in penalised basis coefficients $\theta_{p}$ and Penalised-splines
(P-splines) (Tibshirani, 1996). To obtain P-splines in a Bayesian context a
prior distribution for $\theta_{p}$ and a global smoothing parameter
${\lambda^{G}}$ are introduced. Priors for Basis coefficients $\theta_{p}$ are
given as zero-centred multivariate Normal distribution,
$\boldsymbol{\theta}|\lambda\sim\text{MVN}(0,{\lambda^{G}}^{2}K^{-1})$ with
$\boldsymbol{\theta}$ a vector of $\theta_{p}$ terms (Wood, 2017), and $K$ a
$m$-by-$m$ penalty matrix defined as the squared difference on adjacent basis
coefficients (Eilers & Marx, 1996). Although the zero-centred multivariate-
Normal prior for $\theta_{p}$ is proper, $K$ is rank deficient and the prior
for the basis coefficients is now given as an improper multivariate-Normal
distribution,
$\theta_{p}|\lambda^{G}\propto\text{exp}\left(-\frac{1}{2{\lambda^{G}}^{2}}\boldsymbol{\theta}^{\prime}K\boldsymbol{\theta}\right)$
(Lang & Brezger, 2004). Prior densities for the global shrinkage parameter
have been given as $\lambda^{G}\sim\text{Inverse-Gamma}(1,0.005)$, where a
half-Cauchy, $\lambda^{G}\sim\text{C}^{+}(0,1)$ is preferred to prevent the
prior dominating posterior estimates of global shrinkage (Gelman et al.,
2006). The properties of the additional smoothing parameter is as
$\lambda^{G}\to 0$ the relationship becomes a straight line, and as
$\lambda^{G}\to\infty$ the relationship tends to an un-penalised relationship.

###  3.3 Horseshoe Priors and Local Shrinkage

Although global shrinkage may resolve the issue of over-fitting of smoothed
functions, individual basis coefficients $\theta_{p}$ may become too localised
around zero and become dominated by global shrinkage (Faulkner & Minin, 2018).
Horseshoe priors introduce localised shrinkage parameter
($\lambda_{\theta_{p}}$) by decomposing the variance structure for each basis
coefficient prior into local and global shrinkage (${\lambda^{G}_{\theta}}$)
parameters. Where priors for $\theta_{p}$ and $\lambda_{\theta_{p}}$ are given
as $\theta_{p}|\lambda_{\theta_{p}},\lambda^{G}_{\theta}\sim
N(0,\lambda^{2}_{\theta_{p}}{\lambda^{G}_{\theta}}^{2})$ and
$\lambda_{\theta_{p}}\sim C^{+}(0,1)$ (Carvalho et al., 2010). The
introduction of Horseshoe priors allows for $\theta_{p}$ coefficients to
assume larger values which may otherwise be smoothed out with a single global
shrinkage parameter. Estimating the shrinkage coefficient
$\kappa_{\theta_{p}}$ for each $\lambda_{\theta_{p}}$ is given as

$\kappa_{\theta_{p}}=\frac{1}{1+{\lambda^{G}_{\theta}}^{2}\lambda^{2}_{\theta_{p}}}$

(Carvalho et al., 2010). The global smoothing parameter $\lambda^{G}_{\theta}$
is assumed to be 1 for examples presented in this study, and given
$\lambda_{\theta_{p}}\sim C^{+}(0,1)$ the prior of $\kappa_{\theta_{p}}$ can
be shown to be proportional to a $\text{Beta}(0.5,0.5)$. As
$\kappa_{\theta_{p}}\to 0$ there is no shrinkage of coefficient $\theta_{p}$,
and as $\kappa_{\theta_{p}}\to 1$ there is total shrinkage of coefficient
$\theta_{p}$ to zero (Supplementary Figure 1).

###  3.4 Covariance Structures for Global Shrinkage

With the introduction of local shrinkage priors to $\theta_{p}$, the
underlying prior for $\theta_{p}$ has also changed from a zero-centred
multivariate-Normal distribution for all $m$ terms to $m$ zero-centred Normal
distributions. A series of $m$ zero-centred Normal distributions do provide
adequate local and global shrinkage for smoothed function (Faulkner et al.,
2020). However, posterior estimates of $m$ total basis coefficients are often
correlated and covariate shrinkage of basis coefficients is not identified or
regulated with a series of $m$ univariate Horseshoe priors.

Although a covariance structure for a zero-centred multivariate-Normal prior
on basis coefficients is pre-specified for a global shrinkage parameter
$\lambda^{G}$ given by penalty matrix $K$ (Lang & Brezger, 2004), $K$ is rank
deficient and often specified as the squared difference on adjacent basis
coefficients. An Inverse-Wishart prior on the covariance matrix of the zero-
centred multivariate-Normal prior is a proper prior for global shrinkage with
a covariance structure of basis coefficients (Kass & Natarajan, 2006). Now
$\boldsymbol{\theta}|\Sigma_{\theta},\psi_{\theta}\sim\text{MVN}(0,\Sigma_{\theta})$
and $\Sigma_{\theta}|\psi_{\theta}\sim\text{Inv-
Wishart}(I_{m},\psi_{\theta})$, with $I_{m}$ the identity scale matrix for $m$
coefficients and $\psi_{\theta}$ defined the degrees of freedom. As
$\psi_{\theta}\to\infty$ the Inverse-Wishart prior contributes to more
shrinkage than the Horseshoe prior (Supplementary Figure 1), this results in a
shrinking of the prior sample space for $\theta_{p}$ (Supplementary Figure 2).
Therefore, there is a lack of independence between the number of $\theta_{p}$
and the total shrinkage of basis coefficients, as the degrees of freedom for
the Inverse-Wishart distribution needs to be larger than the number of total
$m$ basis functions. In this study we set $\psi_{\theta}$ to be sufficiently
large at 20 degrees of freedom. $\Sigma_{\theta}$ can be considered as the
generalised covariance form of global shrinkage for basis coefficients without
the need of prior specification of a penalty matrix $K$. With the introduction
of a covariance structure to the shrinkage parameter the generalised covariate
shrinkage coefficient is given as

$\kappa_{\theta_{pq}}=\frac{1}{1+\Sigma^{*}_{\theta_{pq}}}.$

$\Sigma^{*}$ is given as the absolute value of covariance matrix $\Sigma$ as
the prior sampling space for $\theta$ is symmetric. However, $\Sigma_{\theta}$
is not specified with local shrinkage parameters and further considerations of
the decomposition of the covariance structure should be considered to allow
for localised shrinkage and co-shrinkage.

###  3.5 Correlation Structures for Local Shrinkage

For a multivariate Horseshoe prior, consider the vector of basis coefficients
and local shrinkage parameters as $\boldsymbol{\theta}$ and
$\boldsymbol{\lambda_{\theta}}$. The Horseshoe prior for incorporating a
covariance matrix for spline basis coefficients decomposes $\Sigma_{\theta}$
into local and global shrinkage terms along with a correlation structure
$\boldsymbol{\theta}|\boldsymbol{\lambda_{\theta}},\lambda^{G}_{\theta},\Omega_{\theta}\sim\text{MVN}(0,{\lambda^{G}_{\theta}}^{2}\boldsymbol{\lambda_{\theta}}\boldsymbol{\lambda}^{T}_{\boldsymbol{\theta}}\Omega_{\theta})$
(Barnard et al., 2000). With $\Omega_{\theta}$ given as the basis coefficient
correlation matrix with the prior covariance matrix of $\boldsymbol{\theta}$
given as
${\lambda^{G}_{\theta}}^{2}\boldsymbol{\lambda_{\theta}}\boldsymbol{\lambda}^{T}_{\boldsymbol{\theta}}\Omega_{\theta}$.
The prior covariance matrix can be further decomposed into
$\Sigma_{\theta}={\lambda^{G}_{\theta}}^{2}\boldsymbol{\lambda_{\theta}}\boldsymbol{\lambda}^{T}_{\boldsymbol{\theta}}L_{\theta}L^{{}^{\prime}}_{\theta}$.
$L_{\theta}$ is the Cholesky decomposition of the correlation matrix
$\Omega_{\theta}$. The prior for $\Omega_{\theta}$ is given as an
LKJ($\phi_{\theta}$) distribution (Lewandowski et al., 2009) and shape
parameter $\phi_{\theta}>0$. With $\phi_{\theta}=1$ the prior is uniform over
all possible correlation values, $0<\phi_{\theta}<1$ the prior has higher
densities around -1 and 1 with a wider ridge along bivariate parameter
estimates, and as $\phi_{\theta}\to\infty$ the prior is localised around zero
where there is a narrower ridge along bivariate parameter estimates (Figure
2). To regulate $\phi_{\theta}$ estimates, a half-Cauchy$(0,1)$ prior is
assumed. As the correlation between basis functions increases the effect of
covariate shrinkage decreases, whereas as the correlation between coefficients
decreases this introduces more covariate shrinkage (Figure 2). By decomposing
the multivariate-Normal covariance matrix into local shrinkage parameters with
a correlation structure there is more capacity to regulate covariate shrinkage
compared to assuming a Inverse-Wishart prior on the covariance matrix with 20
degrees of freedom (Figure 2 & Supplementary Figure 2).

##  4 Gamma Model

###  4.1 Mean and Variance Parameterisation

Here we specify the reparameterisation of the Gamma model to include TV-AR and
TV-GARCH terms. Consider a non-negative continuous response variable
$Y\sim\text{Gamma}(\alpha,\nu)$ with shape $\alpha$ and rate $\nu$ parameters.
Mean parameterisation of Gamma likelihood is given as
$\mu_{i}=\alpha/\nu_{i}$. With mean estimates of the Gamma likelihood given by
a linear combination of predictors $x_{ij}$ and $j$ model parameters
$\beta_{j}$ via a link function and residual error ($\varepsilon_{i}$), such
that

$\eta_{i}=\sum_{j}x_{ij}\beta_{j}+\varepsilon_{i}.$

The link function for the Gamma distribution is
$\eta_{i}=\sum_{j}x_{ij}\beta_{j}$ such that $\mu_{i}=e^{\eta_{i}}$ to ensure
positive estimates of the response (Dobson, 1990). Estimates of $\nu$ are
given as $\nu_{i}=\alpha/\mu_{i}$ and the likelihood function is
reparameterised by mean estimates $\mu$ and shape parameter $\alpha$, the
likelihood is given as $p(y|\mu,\alpha)$.

The variance from a Gamma process is given as
$\sigma^{2}_{i}=\mu^{2}_{i}/\alpha$ (McCullagh & Nelder, 1983), and $1/\alpha$
is the measure of overall dispersion. Variance structures can be modelled by a
smoothed non-parametric function $g(x_{i}|\tau)$ and basis coefficients for
the variance function $\tau$ (Rice & Silverman, 1991; Silverman, 1985). With
$\sigma^{2}_{i}=g(x_{i}|\tau)$, the shape parameter $\alpha$ is given by
$\alpha_{i}=\sigma^{2}_{i}/\mu^{2}_{i}$ (Addy et al., 2022). The model
likelihood can now be given in terms of mean and variance estimates
$p(y|\mu,\sigma^{2})$.

###  4.2 Time-Varying Autoregression

For an AR(1) process predictions at time $y_{t}$ become
$y_{t}=ay_{t-1}+\varepsilon_{t}$ with predictions based upon the previous
year. TV-AR processes allows for the investigation of how lag-$1$ processes
vary over time. In this manuscript we focus on a TV-AR(1) processes as
grassland herbage yields are know to have a lag-one autoregressive effect. TV-
AR processes can be represented using spline basis functions (Zhang & Wu,
2012) modelling the autoregressive terms as a smoothed non-parametric
relationship through time. Mean estimates from the link function for Gamma
model with regression terms and a TV-AR(1) process becomes

$\eta_{it}=\sum_{j}x_{ij}\beta_{j}+a_{t}y_{t-1}+\varepsilon_{it},$

with $a_{t}=\sum_{p}a_{p}B_{p}(t)$, representing the smooth relationship
between the autoregressive lag-one term and time via TV-AR(1) basis
coefficients $a_{p}$.

###  4.3 Time-Varying GARCH Process

An ARCH process of the $1$st order occurs when the variance of the model
through its squared residual error $\varepsilon_{t}^{2}$ is autoregressive at
lag $1$ (Engle, 1982). A $1$st order GARCH process also considers the moving
average process of the error variance $\sigma^{2}$ along the time-series at
lag $1$ (Bollerslev, 1986). A time-varying GARCH($1$, $1$) process can be
estimated via smoothed non-parametric functions (Karmakar & Roy, 2021) to
identify time-varying processes on the models error variance. In this study we
only consider a TV-GARCH(1, 1) process formally represented as,

$\sigma^{2}_{t}=e^{\tau_{0}+b_{t}\varepsilon^{2}_{t-1}+c_{t}\sigma^{2}_{t-1}}.$

TV-GARCH terms are given as smooth functions, $b_{t}=\sum_{p}b_{p}B_{p}(t)$
and $c_{t}=\sum_{p}c_{p}B_{p}(t)$, and $\tau_{0}$ is the variance intercept.
The GARCH process presented in this study considers heterogeneity in model
residuals conditional on mean estimates.

##  5 Simulation Study

For all simulations within this paper, each non-parametric function was
estimated from 15 basis coefficients from 15 equally spaced Gaussian
dimensions. Standardised residual error at time $t$ was used to model the TV-
ARCH process due to the shifting of the model variance over time, with
$\varepsilon^{s}_{t}=\varepsilon_{t}/\sqrt{\sigma^{2}_{t}}$. A two-step model
fitting process was used, first the TV-AR(1) processes was modelled followed
by the TV-GARCH(1, 1) process. This follows the procedure to model GARCH
processes from the RStan package (Stan Development Team, 2021).

###  5.1 Individual TV-AR(1), TV-ARCH(1) and TV-GARCH(0, 1) Processes

Three individual TV-AR(1), TV-ARCH(1) and TV-GARCH(0, 1) processes were
simulated from 1000 time points with time-varying function

$f(t)=0.45\text{sin}(1.25\times 10^{-2}t)+(5\times 10^{-4}t),$

and $Y\sim\text{Gamma}(\mu,\sigma^{2})$. Although no other time-varying
process was used for each example, we still estimated TV-AR(1) and TV-GARCH(1,
1) processes. This was to assess the stability in the the model approach and
how different prior specifications would estimate a time-varying process even
if one was not present. The fitted model across all prior specifications on
basis functions are given in Supplementary Figures 3-5.

Across all three simulations all three prior specifications for basis
coefficients (Inverse-Wishart prior on a zero-centred multivariate-Normal
covariance matrix, Horseshoe and multivariate Horseshoe) gave good estimates
for the non-parametric time-varying processes. The Inverse-Wishart prior on
the covariance matrix generally gave slightly wigglier estimates, with
Horseshoe and multivariate Horseshoe priors on basis coefficients resulting in
similar smoothed relationships, where the multivariate Horseshoe prior allowed
for regulation of basis coefficients via the LKJ-distribution. Although the
same function was used to simulate all time-varying processes, there was more
variation associated with the time-varying ARCH(1) and GARCH(0, 1) processes
compared to the TV-AR(1) process. Therefore, there is more uncertainty on
variance estimates and the need for better specification on basis function
regulation.

All prior specifications on basis coefficients managed to adequately identify
no relationship for the TV-AR(1) process for the TV-ARCH(1) and TV-GARCH(1)
simulation. However, the use of the Inverse-Wishart prior provided wiggiler
estimates for the TV-ARCH(1) and TV-GARCH(1) when no processes were simulated
compared to the Horseshoe or multivariate Horseshoe prior specifications. This
shows how strict co-regulation can dominate smooth non-parametric
relationships on variance terms and how decomposing the prior covariance
matrix of basis coefficients can adjust for strict co-regulation.

###  5.2 Joint TV-AR(1) and TV-GARCH(1, 1) Process

A TV-AR(1) and TV-GARCH(1, 1) process was simulated for 1000 time points from
a mean-variance parameterised Gamma model $Y\sim\text{Gamma}(\mu,\sigma^{2})$,
with

$\displaystyle a_{t}$
$\displaystyle=\text{sin}\left(\frac{t+10}{75}\right)\left(\frac{50}{t+100}\right),$
$\displaystyle b_{t}$
$\displaystyle=0.5\text{exp}\left(-\left(\frac{x-500}{200}\right)^{10}\right)$
$\displaystyle c_{t}$ $\displaystyle=2.7\times 10^{-9}(t-350)^{3}+5\times
10^{-5}t.$

Where,

$\mu_{t}=e^{\mu_{0}+a_{t}y_{t-1}}$

and

$\sigma^{2}=e^{\tau_{0}+b_{t}\varepsilon^{2}_{t-1}+c_{t}\sigma^{2}_{t-1}},$

with mean and variance intercepts $\mu_{0}=3$ and $\tau_{0}=2.25$. A selection
of different functions for $a_{t}$, $b_{t}$ and $c_{t}$ were selected to
assess the capacity of basis coefficients to deal with over-fitting through
the decomposing of priors into local and global shrinkage terms with a
correlation structure.

The LOOIC (Vehtari et al., 2017) was chosen to compare the posterior
predictive fit of each model. There was no overwhelming differences in the
posterior predictive fit between the models with Inverse-Wishart, Horseshoe or
multivariate Horseshoe (Supplementary Figure 6), suggesting good model fits
across all functions. There were no large differences between the use of
Horseshoe and multivariate Horseshoe priors in the estimation of the TV-
ARCH(1) process for multivariate Horseshoe priors (Figure 3). The use of an
Inverse-Wishart prior on the covariance matrix of basis coefficients tended to
provide wigglier estimates for the TV-ARCH and TV-GARCH functions and a slight
lack of smoothing compared to other prior specifications.

##  6 Park Grass Experiment

###  6.1 Model

Herbage yields from 1966 to 2018 of the Park Grass experiment were modelled
using a mean-variance parameterised Gamma model specified in Section 4 with
$Y\sim\text{Gamma}(\mu,\sigma^{2})$. Mean estimates are given as

$\mu_{jt}=e^{\sum_{j}x_{ij}\beta_{j}+a_{t}y_{t-1}}$

with $\beta_{j}$ a list of seasonally summarised variables presented in Table
1 estimating the effect of climate change on herbage production from 1966 to
2018. The TV-AR(1) process is estimated as a smoothed function $a_{t}$. To
deal with heteroskedastic errors across mean estimates and time-varying
generalised autoregressive conditional heteroskedasticity, the model variance
was estimated as

$\sigma^{2}_{it}=e^{\tau_{i}+b_{t}{\varepsilon^{s}}^{2}_{it-1}+c_{t}\sigma^{2}_{it-1}}.$

The TV-GARCH(1, 1) process was estimated as a smooth function $b_{t}$ and
$c_{t}$. The mean-variance relationship was modelled by a smoothed function
$\tau_{i}=\sum_{p}\tau_{p}B_{p}(\mu_{i})$, with $\tau_{p}$ basis coefficients
and $B_{p}$ basis dimensions with multivariate Horseshoe priors on smooth
terms. It should be noted that the mean-variance relationship and the TV-GARCH
process are modelled from the mean estimate $\mu_{jt}$ given above. The full
list of multivariate Horseshoe priors for the non-parametric mean-variance,
TV-AR(1) and TV-GARCH(1, 1) functions of the Park Grass herbage yield model
are given below.

$\displaystyle\boldsymbol{\tau}|\boldsymbol{\lambda_{\tau}},\lambda_{\tau}^{G},\Omega_{\tau}$
$\displaystyle\sim\text{MVN}(0,{\lambda^{G}_{\tau}}^{2}\boldsymbol{\lambda_{\tau}}\boldsymbol{\lambda}^{T}_{\boldsymbol{\tau}}\Omega_{\tau})$
$\displaystyle\boldsymbol{a}|\boldsymbol{\lambda_{a}},\lambda_{a}^{G},\Omega_{a}$
$\displaystyle\sim\text{MVN}(0,{\lambda^{G}_{a}}^{2}\boldsymbol{\lambda_{a}}\boldsymbol{\lambda}^{T}_{\boldsymbol{a}}\Omega_{a})$
$\displaystyle\boldsymbol{b}|\boldsymbol{\lambda_{b}},\lambda_{b}^{G},\Omega_{b}$
$\displaystyle\sim\text{MVN}(0,{\lambda^{G}_{b}}^{2}\boldsymbol{\lambda_{b}}\boldsymbol{\lambda}^{T}_{\boldsymbol{b}}\Omega_{b})$
$\displaystyle\boldsymbol{c}|\boldsymbol{\lambda_{c}},\lambda_{c}^{G},\Omega_{c}$
$\displaystyle\sim\text{MVN}(0,{\lambda^{G}_{c}}^{2}\boldsymbol{\lambda_{c}}\boldsymbol{\lambda}^{T}_{\boldsymbol{c}}\Omega_{c})$
$\displaystyle\lambda_{\tau_{p}}$ $\displaystyle\sim C^{+}(0,1)$
$\displaystyle\lambda_{a_{p}}$ $\displaystyle\sim C^{+}(0,1)$
$\displaystyle\lambda_{b_{p}}$ $\displaystyle\sim C^{+}(0,1)$
$\displaystyle\lambda_{c_{p}}$ $\displaystyle\sim C^{+}(0,1)$
$\displaystyle\lambda^{G}_{\tau}$ $\displaystyle\sim C^{+}(0,1)$
$\displaystyle\lambda^{G}_{a}$ $\displaystyle\sim C^{+}(0,1)$
$\displaystyle\lambda^{G}_{b}$ $\displaystyle\sim C^{+}(0,1)$
$\displaystyle\lambda^{G}_{c}$ $\displaystyle\sim C^{+}(0,1)$
$\displaystyle\Omega_{\tau}$ $\displaystyle\sim\text{LKJ}(\phi_{\tau})$
$\displaystyle\Omega_{a}$ $\displaystyle\sim\text{LKJ}(\phi_{a})$
$\displaystyle\Omega_{b}$ $\displaystyle\sim\text{LKJ}(\phi_{b})$
$\displaystyle\Omega_{c}$ $\displaystyle\sim\text{LKJ}(\phi_{c})$
$\displaystyle\phi_{\tau}$ $\displaystyle\sim C^{+}(0,1)$
$\displaystyle\phi_{a}$ $\displaystyle\sim C^{+}(0,1)$ $\displaystyle\phi_{b}$
$\displaystyle\sim C^{+}(0,1)$ $\displaystyle\phi_{c}$ $\displaystyle\sim
C^{+}(0,1)$

$\lambda_{\tau}^{G}$, $\lambda_{a_{t}}^{G}$, $\lambda_{b_{t}}^{G}$ and
$\lambda_{c_{t}}^{G}$ are the corresponding global smoothing parameters for
each non-parametric function; $\boldsymbol{\tau}$, $\boldsymbol{a}$,
$\boldsymbol{b}$ and $\boldsymbol{c}$ are a list of estimated terms from a
non-parametric function; and $\Omega_{\tau}$, $\Omega_{a}$, $\Omega_{b}$, and
$\Omega_{c}$ are the estimated correlation matrix sampled from a respective
LKJ(${\phi_{\tau}}$), LKJ($\phi_{a}$), LKJ($\phi_{b}$), and LKJ($\phi_{c}$)
prior. Comparisons of the Inverse-Wishart and Horseshoe prior specifications
for basis coefficients from smooth functions $\tau_{i}$, $a_{t}$, $b_{t}$ and
$c_{t}$ are also given (Figure 4).

###  6.2 Climate Change and Time-Varying Autoregression

From 1966 to 2018, warmer winter temperatures were negatively associated with
lower herbage yields on Park Grass (Table 1). There was a slight positive
linear association with total rainfall in the autumn, winter and summer and an
increase in herbage yields (Table 1). However, there was an inverse polynomial
effect of spring weather (both total rainfall and mean temperature) on herbage
yield. Suggesting too much spring rainfall or warmer spring temperatures would
result in lower herbage yields. The TV-AR(1) process increased in the 1960s,
mid-1980s to late-1990s, and from the late-2000 to 2018, and decreased in the
mid-1970s to 1980s, and mid-1990s to late-2000s (Figure 4). Overall there was
a declining trend of the AR(1) process over time. The estimation of the TV-
AR(1) was consistent across all prior specifications and shows the
autoregressive process a lack of stationarity of the autoregressive process.

###  6.3 Mean-Variance and Time-Varying Variance

The mean-variance relationship was successfully modelled by a non-parametric
smoothed function along with the TV-GARCH(1, 1) process (Figure 4). As herbage
yields increase the variation around herbage yields also increases, therefore
suggesting more yield instability (see yield plots in Figure 1).

The estimated time-varying ARCH(1) process using Horseshoe priors and
multivariate Horseshoe priors declined in the 1980s, and increased from the
1990s to the early 2000s. For the Inverse-Wishart prior specification of the
covariance matrix of basis coefficients the estimated TV-ARCH(1) was wigglier
with more drastic changes (Figure 4). Slight changes in the TV-ARCH(1) process
between 2010 and 2018 were not identified using Horseshoe or multivariate
Horseshoe priors, but were identified using Inverse-Wishart priors on the
covariance matrix of basis coefficients. Given the slight overfitting of the
TV-ARCH(1) process from the first simulation study (Section 5.1), the
regulation from the Horseshoe and multivariate Horseshoe priors of basis
coefficients prevented the slight overfitting of the TV-ARCH(1) term compared
to the Inverse-Wishart prior. The time-varying GARCH process was similar to
the ARCH process with a slight increase in variance between the 1990s and
mid-2000s. This was identified across all prior specifications on basis
coefficients (Figure 4), with Inverse-Wishart prior estimating a slightly
wigglier function. Therefore, with no real differences between the TV-AR(1)
process across prior specifications, the use of Horseshoe and multivariate
Horseshoe priors for basis coefficients was required for preventing slight
overfitting of time-varying variance terms, with no large differences in the
posterior predictive fit across all prior specifications (Supplementary Figure
7).

Model validation of the two-step fitting procedure for the Park Grass data
(TV-AR(1) and then TV-GARCH(1, 1)) follows the Bayesian workflow (Gelman et
al., 2020; Gabry et al., 2019). The Posterior Predictive Distribution (PPD)
for mean estimates from the TV-AR(1) model represented the data well for all
prior specifications (Supplementary Figure 8). Similarity for standard
deviation estimates for the TV-GARCH(1, 1) model, the PPD of the variance
estimates in the Park Grass data were well represented (Supplementary Figure
8). Further diagnostics for LOOIC sampling and Pareto-smoothed importance
sampling (PSIS) are given in Supplementary Figure 9, with most observations
below the Pareto 0.7 threshold (Vehtari et al., 2017).

##  7 Discussion

We proposed the use of the multivariate Horseshoe prior on non-parametric
basis coefficients when modelling time-varying AR(1) and GARCH(1, 1) processes
of the herbage yield data from the Park Grass Experiment, UK. The use of the
multivariate Horseshoe prior decomposes the covariance matrix of a zero-
centred multivariate-Normal prior into a correlation matrix along with local
and global shrinkage terms, with a LKJ-prior on the correlation matrix along
with half-Cauchy(0, 1) priors on local and global shrinkage terms. The
multivariate Horseshoe prior specification on non-parametric basis
coefficients successfully modelled the mean-variance relationship, TV-AR(1)
and TV-GARCH(1, 1) processes of Park Grass herbage yield from 1966 to 2018. We
introduced and assessed the use of the Inverse-Wishart prior (Kass &
Natarajan, 2006) for a non-decomposed covariance matrix on a zero-centred
Normal prior on basis coefficients as a multivariate approach to global
shrinkage of basis coefficients. Compared to previously defined P-splines
(Lang & Brezger, 2004) where a user defined penalty matrix $K$ on the
covariance structure of zero-centred multivariate-Normal distribution is used,
the Inverse-Wishart prior is now proper and not rank deficient. There were no
large differences in the parameterisation of the TV-AR process across the
three different prior specifications (Inverse-Wishart, Horseshoe and
multivariate Horseshoe). However, differences in the prior specifications were
more pronounced in the estimation of the TV-GARCH(1, 1) process suggesting the
co-regulation of basis coefficients from the Inverse-Wishart prior provided
less smoothing of time-varying autocorrelation on variance terms. The
decomposition of a variance matrix was suggested by Barnard et al. (2000) and
in the application to time-varying processes we agree this decomposition
provides more flexibility in covariate shrinkage than the Inverse-Wishart
prior. Further work could be proposed into the use of the LKJ-prior on the
correlation matrix of a regulated Horseshoe priors (Piironen & Vehtari, 2017),
and how co-regulation could be achieved. Other prior distributions may also be
used for the basis coefficients correlation terms such as the Beta
distribution (Gokhale & Press, 1982). However, this would require
$\frac{3m!}{2!(m-2)}$ priors and with two shape parameters each which could
prove computationally intensive.

For the Park Grass data from 1966 to 2018 we used a generalised mean-variance
parameterisation of the Gamma function to model the TV-AR(1) and TV-GARCH(1,
1) processes. TV-AR(1) and TV-GARCH(1, 1) processes have been used previously
for count GLMs (Roy & Karmakar, 2021), however we generalised this process for
continuous non-negative data via a Gamma model. From the analysis of the Park
Grass herbage yield data we conclude similar results to (Addy et al., 2021;
Cashen, 1947) that warmer temperatures in spring and winter resulted in lower
herbage yields. In both Addy et al. (2022) and this study, we observed a
consistently positive autocorrelation between yields, indicating that yield-
related ecological processes such as seed production and/or energy storage
play a role in determining yields in the following year (Ma et al., 2021).
However, in this study, by modelling a TV-AR(1) process of herbage yield we
also found the autoregressive process to be non-stationary, but always
positive, alongside the influence of climate change and warmer temperatures
(Figure 4). Although the mechanism behind the time-varying pattern in yield
autocorrelation is not known, it is of interest especially for forecasting
future herbage yields into the mid and late-21st century (Addy et al., 2022).

By formalising the modelling procedure to include time-varying stochastic time
series processes on Park Grass, the TV-GARCH(1, 1) process showed how the
variance around the model prediction of the Park Grass herbage yield data are
autoregressive and have slightly increased between the 1990s and early-2000s,
confirming previous findings in yield stability changing during the 1990s
(Macholdt et al., 2023; Dodd et al., 1994). With the mean-variance
relationship, we also found that changes in yield stability over time were
more severe for higher-yield sections of the Park Grass Experiment. Macholdt
et al. (2023) observed a peak in yield variance in the 1990s, at least
partially explained by an increase in summer water stress and temperatures
during this decade. In contrast, the identification of the TV-GARCH process
suggests that variance increased in the 1990s and remained high from the 2000s
onwards (Figure 4), and meteorological trends from the Park Grass site do
suggest a sustained increase in temperatures from the 1990s onwards (Addy et
al., 2022). Zhang et al. (2018) observed that higher temperatures could
decrease grassland species richness and consequently yield stability in
northern China, so higher temperatures may have had a similar effect on Park
Grass. Another ecologically significant event that occurred during the 1990s
was a drop in nitrogen and sulphur deposition, which was associated with an
increase in soil pH and an increased proportion of legume species in the
herbage cuts (Storkey et al., 2015; Macholdt et al., 2023), and this may also
have affected yield stability.

Our conclusions differ slightly to previous work on identifying the drivers of
yield stability in Park Grass, with Macholdt et al. (2023) observing an
increase in variance in the 1990s but a subsequent decrease, whilst we
identified a sustained increase in variance from the 1990s onwards. We
modelled different treatments (only two were consistent), and we also modelled
yield variance after estimating mean predictions from weather variables and a
TV-AR(1) process. Our finding show that Park Grass yield stability has
followed a consistently different pattern since the 1990s and this has
important implications for understanding how grasslands respond to climate
change. However, it is important to note that the range of values from the
time-varying GARCH(1, 1) process was smaller compared to the mean-variance
process (Figure 4). Overall, our paper demonstrates that variance functions
are highly sensitive and dependent on regulation, and considerations should be
given when trying to attribute changes in variance over time to other
stochastic processes such as weather.

## References

  * (1)
  * Addy et al. (2021) Addy, J. W., Ellis, R. H., Macdonald, A. J., Semenov, M. A. & Mead, A. (2021), ‘Changes in agricultural climate in south-eastern england from 1892 to 2016 and differences in cereal and permanent grassland yield’, Agricultural and Forest Meteorology 308, 108560.
  * Addy et al. (2022) Addy, J. W., Ellis, R. H., MacLaren, C., Macdonald, A. J., Semenov, M. A. & Mead, A. (2022), ‘A heteroskedastic model of Park Grass spring hay yields in response to weather suggests continuing yield decline with climate change in future decades’, Journal of the Royal Society Interface 19(193).
  * Barnard et al. (2000) Barnard, J., McCulloch, R. & Meng, X.-L. (2000), ‘Modeling covariance matrices in terms of standard deviations and correlations, with application to shrinkage’, Statistica Sinica pp. 1281–1311.
  * Bengtsson et al. (2019) Bengtsson, J., Bullock, J., Egoh, B., Everson, C., Everson, T., O’connor, T., O’farrell, P., Smith, H. & Lindborg, R. (2019), ‘Grasslands—more important for ecosystem services than you might think’, Ecosphere 10(2), e02582.
  * Bollerslev (1986) Bollerslev, T. (1986), ‘Generalized autoregressive conditional heteroskedasticity’, Journal of Econometrics 31(3), 307–327.
  * Bova et al. (2021) Bova, S., Rosenthal, Y., Liu, Z., Godad, S. P. & Yan, M. (2021), ‘Seasonal origin of the thermal maxima at the holocene and the last interglacial’, Nature 589(7843), 548–553.
  * Carvalho et al. (2010) Carvalho, C. M., Polson, N. G. & Scott, J. G. (2010), ‘The horseshoe estimator for sparse signals’, Biometrika 97(2), 465–480.
  * Cashen (1947) Cashen, R. O. (1947), ‘The influence of rainfall on the yield and botanical composition of permanent grass at Rothamsted’, The Journal of Agricultural Science 37(1), 1–10.
  * Dobson (1990) Dobson, A. J. (1990), An Introduction to Generalized Linear Models, CRC press.
  * Dodd et al. (1994) Dodd, M. E., Silvertown, J., McConway, K., Potts, J. & Crawley, M. (1994), ‘Stability in the plant communities of the park grass experiment: the relationships between species richness, soil pH and biomass variability’, Philosophical Transactions of the Royal Society of London. Series B: Biological Sciences 346(1316), 185–193.
  * Eilers & Marx (1996) Eilers, P. H. & Marx, B. D. (1996), ‘Flexible smoothing with B-splines and penalties’, Statistical Science 11(2), 89–121.
  * Engle (1982) Engle, R. F. (1982), ‘Autoregressive conditional heteroscedasticity with estimates of the variance of united kingdom inflation’, Econometrica pp. 987–1007.
  * Faulkner et al. (2020) Faulkner, J. R., Magee, A. F., Shapiro, B. & Minin, V. N. (2020), ‘Horseshoe-based bayesian nonparametric estimation of effective population size trajectories’, Biometrics 76(3), 677–690.
  * Faulkner & Minin (2018) Faulkner, J. R. & Minin, V. N. (2018), ‘Locally adaptive smoothing with markov random fields and shrinkage priors’, Bayesian Analysis 13(1), 225.
  * Gabry et al. (2019) Gabry, J., Simpson, D., Vehtari, A., Betancourt, M. & Gelman, A. (2019), ‘Visualization in Bayesian workflow’, Journal of the Royal Statistical Society Series A: Statistics in Society 182(2), 389–402.
  * Gelman et al. (2020) Gelman, A., Vehtari, A., Simpson, D., Margossian, C. C., Carpenter, B., Yao, Y., Kennedy, L., Gabry, J., Bürkner, P.-C. & Modrák, M. (2020), ‘Bayesian workflow’, arXiv preprint arXiv:2011.01808 .
  * Gelman et al. (2006) Gelman, A. et al. (2006), ‘Prior distributions for variance parameters in hierarchical models (comment on article by Browne and Draper)’, Bayesian Analysis 1(3), 515–534.
  * Giudici et al. (2023) Giudici, P., Tarantino, B. & Roy, A. (2023), ‘Bayesian time-varying autoregressive models of covid-19 epidemics’, Biometrical Journal 65(1), 2200054.
  * Gokhale & Press (1982) Gokhale, D. & Press, S. J. (1982), ‘Assessment of a prior distribution for the correlation coefficient in a bivariate Normal distribution’, Journal of the Royal Statistical Society: Series A (General) 145(2), 237–249.
  * Hoover et al. (1998) Hoover, D. R., Rice, J. A., Wu, C. O. & Yang, L.-P. (1998), ‘Nonparametric smoothing estimates of time-varying coefficient models with longitudinal data’, Biometrika 85(4), 809–822.
  * IPCC (2023) IPCC (2023), ‘Synthesis Report of the IPCC Sixth Assessment Report’.
https://www.ipcc.ch/report/ar6/syr/

  * Jenkinson et al. (1994) Jenkinson, D., Potts, J., Perry, J., Barnett, V., Coleman, K. & Johnston, A. (1994), ‘Trends in herbage yields over the last century on the Rothamsted long-term continuous hay experiment’, The Journal of Agricultural Science 122(3), 365–374.
  * Karmakar & Roy (2021) Karmakar, S. & Roy, A. (2021), ‘Bayesian modelling of time-varying conditional heteroscedasticity’, Bayesian Analysis 16(4), 1157–1185.
  * Kass & Natarajan (2006) Kass, R. E. & Natarajan, R. (2006), ‘A default conjugate prior for variance components in generalized linear mixed models (comment on article by Browne and Draper)’, Bayesian Analysis 1(3), 535–542.
  * Lang & Brezger (2004) Lang, S. & Brezger, A. (2004), ‘Bayesian P-splines’, Journal of Computational and Graphical Statistics 13(1), 183–212.
  * Lewandowski et al. (2009) Lewandowski, D., Kurowicka, D. & Joe, H. (2009), ‘Generating random correlation matrices based on vines and extended onion method’, Journal of Multivariate Analysis 100(9), 1989–2001.
  * Little & Rubin (1986) Little, R. J. & Rubin, D. B. (1986), Statistical analysis with missing data, John Wiley & Sons.
  * Ma et al. (2021) Ma, M., Collins, S. L., Ratajczak, Z. & Du, G. (2021), ‘Soil seed banks, alternative stable state theory, and ecosystem resilience’, BioScience 71(7), 697–707.
  * Macdonald (2018) Macdonald, A. (2018), Guide to the Classical and Other Long-Term Experiments, Datasets and Sample Archive, Rothamsted Research.
  * Macholdt et al. (2023) Macholdt, J., Hadasch, S., Macdonald, A., Perryman, S., Piepho, H.-P., Scott, T., Styczen, M. E. & Storkey, J. (2023), ‘Long-term trends in yield variance of temperate managed grassland’, Agronomy for Sustainable Development 43(3), 37.
  * McCullagh & Nelder (1983) McCullagh, P. & Nelder, J, A. (1983), Generalized Linear Models, Chapman and Hall.
  * Piepho (1998) Piepho, H.-P. (1998), ‘Methods for comparing the yield stability of cropping systems’, Journal of Agronomy and Crop Science 180(4), 193–213.
  * Piironen & Vehtari (2017) Piironen, J. & Vehtari, A. (2017), ‘Sparsity information and regularization in the Horseshoe and other shrinkage priors’, Electronic Journal of Statistics 11, 5018–5051.
  * Rice & Silverman (1991) Rice, J. A. & Silverman, B. W. (1991), ‘Estimating the mean and covariance structure nonparametrically when the data are curves’, Journal of the Royal Statistical Society: Series B (Methodological) 53(1), 233–243.
  * Roy & Karmakar (2021) Roy, A. & Karmakar, S. (2021), ‘Time-varying auto-regressive models for count time-series’, Electronic Journal of Statistics 15(1), 2905–2938.
  * Silverman (1985) Silverman, B. W. (1985), ‘Some aspects of the spline smoothing approach to non-parametric regression curve fitting’, Journal of the Royal Statistical Society: Series B (Methodological) 47(1), 1–21.
  * Stan Development Team (2021) Stan Development Team (2021), ‘RStan: the R interface to Stan’.  R package version 2.21.3.
https://mc-stan.org/

  * Storkey et al. (2015) Storkey, J., Macdonald, A., Poulton, P., Scott, T., Köhler, I., Schnyder, H., Goulding, K. & Crawley, M. (2015), ‘Grassland biodiversity bounces back from long-term nitrogen addition’, Nature 528(7582), 401–404.
  * Tibshirani (1996) Tibshirani, R. (1996), ‘Regression shrinkage and selection via the lasso’, Journal of the Royal Statistical Society: Series B (Methodological) 58(1), 267–288.
  * Vehtari et al. (2017) Vehtari, A., Gelman, A. & Gabry, J. (2017), ‘Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC’, Statistics and Computing 27(5), 1413–1432.
  * Willcock et al. (2023) Willcock, S., Cooper, G. S., Addy, J. & Dearing, J. A. (2023), ‘Earlier collapse of anthropocene ecosystems driven by multiple faster and noisier drivers’, Nature Sustainability pp. 1–12.
  * Wood (2017) Wood, S. N. (2017), Generalized Additive Models: an Introduction with R, CRC press.
  * Zhang et al. (2022) Zhang, F., Bennett, J. A., Zhang, B., Zhao, M. & Han, G. (2022), ‘Intensification of disturbance destabilizes productivity through effects on dominant species’, Ecological Indicators 143, 109383.
  * Zhang & Wu (2012) Zhang, T. & Wu, W. B. (2012), ‘Inference of time-varying regression models’, The Annals of Statistics 40(3), 1376–1402.
  * Zhang et al. (2018) Zhang, Y., Loreau, M., He, N., Wang, J., Pan, Q., Bai, Y. & Han, X. (2018), ‘Climate variability decreases species richness and community stability in a temperate grassland’, Oecologia 188, 183–192.

## Figures and Tables

Figure 1: Herbage yields (t dm ha-1) from all Park Grass limed sections plots
14.2, 16, 7, 2.2 and 3 from 1966 to 2018. Limed sections are staggered from no
liming (light grey) to pH maintained at 7 (black).

Figure 2: The distribution of correlation terms from 10000 simulations of the
LKJ-distribution for different $\phi$ values (top row). The bivariate
distribution of basis coefficients from a multivariate Horseshoe distribution
with correlation terms sampled from a prior LKJ-distribution for different
$\phi$ values (middle row). The influence of covariate shrinkage between basis
coefficients with varying correlations sampled from 10000 simulations (bottom
row).

Figure 3: The simulated TV-AR(1) and TV-GARCH(1, 1) processes for a mean-
variance parameterised Gamma model (top row). The fitted TV-AR(1) (second row)
and TV-GARCH(1, 1) (rows three and four) with 80% credible intervals for each
type of basis coefficient prior; Inverse-Wishart($20$) prior on covariance
matrix (green), Horseshoe (blue) and multivariate Horseshoe (red). All
simulated data and functions are given in grey.

Figure 4: The mean-variance function ($\tau_{i}$), TV-AR(1) and TV-GARCH(1, 1)
processes of the mean-variance parameterised Gamma model for the Park Grass
data with 80% credible intervals for each type of basis coefficient prior;
Inverse-Wishart($20$) prior on covariance matrix (green), Horseshoe (blue) and
multivariate Horseshoe (red).

Figure 5: Variance estimates from mean-variance and TV-GARCH(1, 1) process
with multivariate Horseshoe priors of Park Grass limed sections plots 14.2,
16, 7, 2.2 and 3 from 1966 to 2018. Limed sections are staggered from no
liming (light grey) to pH maintained at 7 (black).

Table 1: Weather parameter estimates with standard deviation (SD), 80% credible intervals from the Park Grass regression model. TR and MT refer to Total Rainfall and Mean Temperature with subscript 1 and 2 correspond to linear and polynomial terms. Model Terms | Mean | SD | 10% | 90%
---|---|---|---|---
$\beta_{\text{TR Aut (mm)}}$ | 3.40$\times 10^{-4}$ | 1.120$\times 10^{-4}$ | 1.97$\times 10^{-4}$ | 4.84$\times 10^{-4}$
$\beta_{\text{TR Win (mm)}}$ | 2.56$\times 10^{-4}$ | 1.124$\times 10^{-4}$ | 1.13$\times 10^{-4}$ | 3.98$\times 10^{-4}$
$\beta_{\text{TR Spr}_{1}\text{(mm)}}$ | 1.99 | 0.266 | 1.64 | 2.33
$\beta_{\text{TR Spr}_{2}\text{(mm)}}$ | -1.74 | 0.308 | -2.12 | -1.34
$\beta_{\text{TR Sum (mm)}}$ | 1.74$\times 10^{-3}$ | 0.186$\times 10^{-3}$ | 1.50$\times 10^{-3}$ | 1.99$\times 10^{-3}$
$\beta_{\text{MT Aut (deg C)}}$ | 2.23$\times 10^{-2}$ | 1.056$\times 10^{-2}$ | 0.84$\times 10^{-2}$ | 3.58$\times 10^{-2}$
$\beta_{\text{MT Win (deg C)}}$ | -3.83$\times 10^{-2}$ | 0.815$\times 10^{-2}$ | -4.86$\times 10^{-2}$ | -2.78$\times 10^{-2}$
$\beta_{\text{MTSpr}_{1}\text{(deg C)}}$ | 1.19 | 0.340 | 0.72 | 1.65
$\beta_{\text{MTSpr}_{2}\text{(deg C)}}$ | -1.89 | 0.259 | -2.22 | -1.57
$\beta_{\text{MT Sum (deg C)}}$ | 2.96$\times 10^{-2}$ | 1.162$\times 10^{-2}$ | 1.45$\times 10^{-2}$ | 4.49$\times 10^{-2}$

## Supplementary Materials

Supplementary Figure 1: 10000 samples of the shrinkage coefficient from an
Inverse-Wishart($\psi=20$) and Horseshoe prior. Supplementary Figure 2: 10000
samples of a bivariate Inverse-Wishart prior with $\psi=3,10$ and $20$ (top
row). The covariate shrinkage of 10000 samples from an Inverse-Wishart prior
with $\psi=3,10$ and $20$ (bottom row).

Supplementary Figure 3: The simulated TV-AR(1) processes from simulation
example 5.1 for a mean-variance parameterised Gamma model (top row). The
fitted TV-AR(1) (second row) and TV-GARCH(1, 1) (rows three and four) with 80%
credible intervals for each type of basis coefficient prior; Inverse-
Wishart($20$) prior on covariance matrix (green), Horseshoe (blue) and
multivariate Horseshoe (red). All simulated data and functions are given in
grey.

Supplementary Figure 4: The simulated TV-ARCH(1) process from simulation
example 5.1 for a mean-variance parameterised Gamma model (top row). The
fitted TV-AR(1) (second row) and TV-GARCH(1, 1) (rows three and four) with 80%
credible intervals for each type of basis coefficient prior; Inverse-
Wishart($20$) prior on covariance matrix (green), Horseshoe (blue) and
multivariate Horseshoe (red). All simulated data and functions are given in
grey.

Supplementary Figure 5: The simulated TV-GARCH(1) process from simulation
example 5.1 for a mean-variance parameterised Gamma model (top row). The
fitted TV-AR(1) (second row) and TV-GARCH(1, 1) (rows three and four) with 80%
credible intervals for each type of basis coefficient prior; Inverse-
Wishart($20$) prior on covariance matrix (green), Horseshoe (blue) and
multivariate Horseshoe (red). All simulated data and functions are given in
grey.

Supplementary Figure 6: The LOOIC ($+-$standard error) for the simulation
study for each prior selection on basis coefficients. Supplementary Figure 7:
The LOOIC ($+-$standard error) for the Park Grass example for each prior
selection on basis coefficients.

Supplementary Figure 8: The Posterior Predictive Distribution (PPD) of the
mean estimates of the TV-AR(1) model for Park Grass for all prior
specifications (left). PPD of the standard deviation estimates for the TV-
GARCH(1, 1) model for Park Grass for all prior specifications (right). Both
mean and SD’s were derived from the log-scale of the data.

Supplementary Figure 9: The Pareto-smoothed importance sampling (PSIS) output
from the TV-AR(1) and TV-GARCH(1, 1) models of the Park Grass dataset across
all prior specifications.