Towards Black-Box Parameter Estimation

Amanda Lenzi111 School of Mathematics, University of Edinburgh, Edinburgh EH9
3FD, Scotland, United Kingdom. and Haavard Rue222 Statistics Program, King
Abdullah University of Science and Technology, Thuwal 23955-6900, Saudi
Arabia.

Abstract

Deep learning algorithms have recently been shown to be a successful tool in
estimating parameters of statistical models for which simulation is easy, but
likelihood computation is challenging. This is achieved by sampling a large
number of parameter values from a distribution, which is typically chosen to
be non-informative and cover as much of the parameter space as possible.
However, for high-dimensional and large parameter spaces, covering all
possible reasonable parameter values is infeasible. We propose a new
sequential training procedure that reduces simulation cost and guides
simulations toward the region of high parameter density based on estimates of
the neural network and the observed data. Our following proposal aims to fit
time series models to newly collected data at no cost using a pre-trained
neural network with simulated time series of a fixed length. These approaches
can successfully estimate and quantify the uncertainty of parameters from non-
Gaussian models with complex spatial and temporal dependencies. The success of
our methods is a first step towards a fully flexible automatic black-box
estimation framework.

Keywords: Deep neural networks, intractable likelihoods, sequential, time-
series, simulation

Short title: Black-box Estimation

##  1 Introduction

Statistical modeling consists of first devising stochastic models for
phenomena we want to learn about and then to relate those models to data.
These stochastic models have unknown parameters and the second step boils down
to estimating these parameters from data through the likelihood function.
However, there might be a discrepancy between these two steps, as models for
describing mechanisms aim for scientific adequacy rather than computational
tractability. Indeed, as soon as we move away from Gaussian processes as the
canonical model for dependent data, likelihood computation becomes effectively
impossible, and inference is too complicated for traditional estimation
methods. Consider, for instance, datasets from finance or climate science,
where skewness and jumps are commonly present and calculating the likelihood
in closed form is often impossible, ruling out any numerical likelihood
maximization and Bayesian methods. Yet it is computationally inexpensive to
simulate from those models given parameter values, and the question becomes
whether the simulations look like the data.

Much effort has been directed toward the development of approximate parameter
estimation methods, often referred to as indirect inference (Gourieroux et
al.,, 1993), likelihood-free inference (Grelaud et al.,, 2009; Gutmann and
Corander,, 2016), simulation-based inference (Nickl and Pötscher,, 2010) or
synthetic likelihood (Wood,, 2010); for an overview, see, for example, the
review by Hartig et al., (2011); Cranmer et al., (2020). The typical
assumption by the different methods is that exact likelihood evaluation is
hard to obtain but it is easy to simulate from the model given the parameter
values, and the basic idea is to identify the model parameters which yield
simulated data that resemble the observed data. The most common in this
umbrella is arguably approximate Bayesian computation (ABC) (Fearnhead and
Prangle,, 2012; Frazier et al.,, 2018; Sisson et al.,, 2018), which avoids
evaluating intractable likelihoods by matching summary statistics from
observations with those computed from simulated data based on parameters drawn
from a predefined prior distribution. The likelihood is approximated by the
probability that the condition
$\gamma(x_{\mbox{sim}},x_{\mbox{obs}})<\epsilon$ is satisfied, where $\gamma$
is some distance measure and the value of $\epsilon$ is a trade-off between
sample efficiency and inference quality. In simpler cases, sufficient
statistics are used as they provide all the information in the data, however,
for complex models they are unlike to exist and it is not obvious which
statistics will be most informative. Several works have proposed procedures
for designing summary statistics (Fearnhead and Prangle,, 2012; Jiang et al.,,
2017), and a comparison of likelihood-free methods with and without summary
statistics have been empirically tested in Drovandi and Frazier, (2022).
Despite its popularity, ABC is not scalable to large numbers of observations
since inference for new data requires repeating most steps of the procedure.

A recent line of research on likelihood-free inference uses deep learning to
estimate parameters of statistical models. The first work to propose neural
networks-based estimators in the statistical community focused on parameters
of spatial covariance functions in Gaussian processes (Gerber and Nychka,,
2020). They showed that convolutional neural networks (CNNs) can learn the
mapping between data and parameters and had similar estimation accuracy and a
considerable reduction in computational time compared to classical maximum
likelihood estimators. Based on the same idea, Lenzi et al., (2021) estimated
parameters of models for spatial extremes for which the likelihoods are
intractable and, therefore, MLEs are unavailable. A modified parametric
bootstrap approach was introduced to quantify the uncertainty in these
estimators. Variants to those methods to incorporate replicated data in the
estimation (Sainsbury-Dale et al.,, 2022), irregular spatial data (Sainsbury-
Dale et al.,, 2023), and censoring information (Richards et al.,, 2023) were
successfully recently introduced. Estimation in high dimensions and large
parameter spaces is still an open question since constructing training data
that covers all possible reasonable parameter values becomes quickly
infeasible. An inevitable drawback of the current approaches is the
estimator’s bias towards the parameter region of the training data. Whereas
previous work simulated training data either with parameter values around the
truth or based on inexact likelihood estimates, these methods are expensive
and unrealistic. In this work, we solve the bias issue of previous methods
with an automatic iterative approach that modifies the training data using
arbitrary, dynamically updated distribution parameters until it reaches the
parameter region corresponding to the actual data. The proposed inference
mechanism automatically performs parameter estimation without restrictive
assumptions about the generating process, knowledge from experts, or
computationally expensive preliminary steps.

Various approaches have been proposed for guiding simulations by making use of
the observed data. The sequential training procedure proposed here falls into
a broader class of methods that seek to update the prior distribution to
alleviate the curse of dimensionality that comes with having to explore a
large volume of data when the number of parameters is large. Within ABC
implementations, a Sequential Monte Carlo procedure was used to guide
simulations based on previously accepted parameters (Sisson et al.,, 2007;
Beaumont et al.,, 2009; Bonassi and West,, 2015). Järvenpää et al., (2019);
Lueckmann et al., (2019) updated simulations that reduce the Bayesian
uncertainty in the posterior estimate. A sequential estimation procedure
selected future simulations by proposing parameters from preliminary
approximations to the posterior (Lueckmann et al.,, 2017; Papamakarios and
Murray,, 2016). Sequential approaches have been considered in the context of
neural networks for learning a model of the likelihood in the region of high
posterior density Papamakarios et al., (2019), but only for toy examples and
small datasets. Our approach stands out from previous methods as it is aimed
at point parameter estimation in intractable statistical models with large and
complex datasets. We leverage our statistical knowledge about scaling data and
parameters to reduce computational costs and improve convergence when training
the neural networks. We first illustrate the new approach on a Gaussian toy
example, where the intended coverage probability of the estimator is known,
and further for modeling spatial extremes. Similarly to Papamakarios et al.,
(2019), which trains autoregressive flows on all simulations obtained up to
each round rather than training only with simulations from the latest round,
we broaden the range of the training by reusing training data from previous
iterations. Consequently, most of the training data in later rounds still come
from the most probable parameter regions, and the added simulations from
previous rounds broaden the range of searches at no computational cost. We
leverage the ideas of improving the estimator by updating the prior in
multiple rounds. Since our approach is designed to yield precise point
estimates instead of approximating posterior distributions, the loss function
in the neural network remains independent of the current proposed prior. We
apply the modified parametric bootstrap technique in Lenzi et al., (2021) such
that in each round, the new approach is able to quantify both the bias and
variance of the estimators. Our experimental results show that the sequential
approach dramatically reduces the bias of an initial guess and eventually
approximates the actual parameter quite accurately, even when the initial
training does not contain the truth.

Our next contribution is designed to handle estimation for different data
sizes originating from the same model. This is achieved by training a neural
network only once on an extensive database and replicating the observations as
needed to achieve the training data length before processing them on the pre-
trained neural network. We adjust for the underestimation in the uncertainty
of the replicated series by rescaling its sampling variance to match the
original series variance. We apply this strategy to estimate parameters of
non-Gaussian stochastic volatility models, which are widely used in finance.
In such applications, where early access to results may carry a premium, our
deep neural network (DNN) estimator is particularly advantageous since the
network is trained beforehand, and estimates are obtained instantaneously when
new data becomes available. This example shows that our estimator is well
calibrated, with uncertainty quantification closely matching those from the
state-of-the-art Integrated Nested Laplace Approximation (INLA) approach (Rue
et al.,, 2009).

The remainder of this paper is organized as follows. First, Section 2 outlines
the methodologies we develop for designing training data to train DNNs, along
with some practical considerations. In Section 3, we introduce the
construction of the automatic iterative approach and conduct simulation
studies for Gaussian i.i.d and spatial extremes model, whereas in Section 4,
we describe our unified database approach applied to time series data from an
intractable model. In Section 5, we conclude and summarize avenues for future
research.

##  2 Parameter estimation with DNNs

###  2.1 Background

Consider a dataset $\mathbf{x}_{0}\in\mathbb{R}^{J}$ of observations generated
from $\mathbb{P}\in\mathcal{P}(\mathcal{X})$, where $\mathbb{P}$ is a Lesbegue
measure and $\mathcal{P}(\mathcal{X})$ denote the set of all Borel probability
measures on the sample space $\mathcal{X}$. To describe such a process, it is
common practice to assume a statistical model
$\mathbb{P}_{\bm{\theta}_{0}}\in\mathcal{P}(\mathcal{X})$ with probability
density function $p(\cdot;\bm{\theta}_{0})$ parameterized by a finite number
of parameters $\bm{\theta}_{0}\subset\Theta\in\mathbb{R}^{P}$, which is
estimated using the observations through the log-likelihood function
$l(\bm{\theta}_{0};\mathbf{x}_{0})\equiv\mbox{log}\\{p(\mathbf{x}_{0};\bm{\theta}_{0})\\}$.

Highly structured data coming from a high-dimensional $\mathcal{X}$ are often
related to intractable or computationally demanding likelihoods, but
simulating data from $p$ for given parameters is usually trivial. Recently,
parameter estimation using DNNs have opened doors to solving previously
intractable statistical estimation problems. The key to efficiency is to avoid
altogether learning likelihood functions and directly learn the mapping
between data and parameters through DNNs by carrying out simulations. To
formulate the problem, let $\mathbf{x}_{n}\in\mathbb{R}^{J}$ be a simulated
sample from $p(\cdot;\bm{\theta}_{n})$ with given parameters
$\bm{\theta}_{n}\in\mathbb{R}^{P}$. Then, the mapping from
$\mathbf{x}=(\mathbf{x}_{1},\ldots,\mathbf{x}_{N})^{\top}$ onto
$\bm{\theta}=(\bm{\theta}_{1},\ldots,\bm{\theta}_{N})^{\top}$ is learned by
adjusting the weights w and biases b, denoted by
$\phi=(\textbf{w},\textbf{b})^{\top}$ of a DNN $\mathcal{F}_{\phi}$, such that

$\mathcal{F}_{\phi}:\mathbf{x}\mapsto\bm{\theta};\quad\hat{\bm{\theta}}=\mbox{argmin}_{\phi}\mathcal{L}\\{\bm{\theta},\mathcal{F}_{\phi}(\mathbf{x})\\}.$
(1)

Optimizing (1) with respect to $\phi$ requires the minimization of the loss
function $\mathcal{L}$, which is chosen to reduce the error in prediction for
a given output $\bm{\theta}$ and simulated data $\mathbf{x}$. A popular choice
in regression problems is the mean squared error (MSE)

$\mathcal{L}(\phi;\bm{\theta},\mathbf{x})=\mathbb{E}\\{\bm{\theta}-\mathcal{F}_{\phi}(\mathbf{x})\\}^{2}.$

Often, no closed-form solutions can be derived for the optimization in (1),
and advanced numerical optimizers built around batch gradient descent methods
are employed (Kingma and Ba,, 2014). Finally, once the DNN has been trained,
one can use the estimated $\hat{\phi}$ to plug in $\mathbf{x}_{0}$ into the
trained DNN and retrieve $\mathcal{F}_{\hat{\phi}}(\mathbf{x}_{0})$, which
will then output parameter estimates of interest $\hat{\bm{\theta}}_{0}$.

###  2.2 Transformations to data and parameters

Here, we detail our rationale for choosing transformations to data and
parameters and make the problem more palatable for the DNN. The key is to use
our statistical knowledge of intrinsic data and parameter properties to
leverage estimation. For instance, the quadratic loss in (2.1) is optimal for
outputs with constant mean and variance that are a real-valued function of the
inputs with Gaussian distributed noises. If these assumptions are met, the
estimator $\mathcal{F}$ retains the desired properties, such as the minimum
variance and fast convergence, as the gradient reduces gradually for
relatively small errors (Friedman et al.,, 2001). Therefore, reparametrization
related to scaling should aim for a constant variance for different parameter
values.

When minimizing the loss in (1), $\phi$ is usually initialized to random
values and updated via an optimization algorithm, such as stochastic gradient
descent based on training data. As is usually the case, the geometry of the
surface that has to be optimized will be complicated and smooth due to an
ample search space and noisy data, and using the raw data will likely result
in slow and unstable convergence. To remedy this problem and improve the
algorithm’s stability, one should aim for properties such as symmetric and
unbounded distributions, orthogonal parameters, and constant Fisher
information. The logarithm and square root transformations often used in time
series problems are examples of desirable change that affects the distribution
shape by reducing skewness, stabilizing the variance, and simultaneously
avoiding boundaries. Parameterization with meaningful interpretations such as
mean and variance should be preferred over directly using distribution
parameters.

In Sections 3 and 4, we will use transformation within a DNN pipeline to
estimate parameters of models for Gaussian data, spatial extremes, and non-
Gaussian stochastic volatility models. We show that whereas these precautions
are helpful in simple Gaussian examples, they are indispensable in complex
models such as for spatial extremes.

###  2.3 Designing training data

Recall that the first step for optimizing (1) is generating pairs of training
data $(\bm{\theta}_{n},\mathbf{x}_{n})^{N}_{n=1}$. The main challenge here is
to generate training data that correspond to configurations covering the
parameter domain of the observations, which is unknown. Since $\bm{\Theta}$ is
often unbounded and it is impossible to simulate over the entire domain. Not
introducing an appropriate structure or prior scientific knowledge will lead
to inaccurate training data and, thus, erroneous estimation.

In the context of intractable likelihoods for which MLEs are unavailable,
Lenzi et al., (2021) proposed to simulate training data based on informative
parameter estimates from approximate maximum likelihood methods fit to spatial
extremes data. However, obtaining these estimates for every new dataset
becomes problematic if likelihood estimation is slow or not feasible in the
first place. Here, we propose two different strategies to deal with this
challenge. The intuition behind these method goes as follow:

  1. (A)

A fully automatic iterative approach: Promising regions in $\bm{\Theta}$ are
found sequentially. The DNN is initially trained with $\bm{\theta}$ based on a
crude guess (e.g., from a simpler model). The trained DNN then receives
$\mathbf{x}_{0}$ and outputs $\hat{\bm{\theta}}_{0}$, which is used to
simulate bootstrapping samples $\mathbf{x}_{b}\sim p(;\hat{\bm{\theta}}_{0})$.
Next, $\mathbf{x}_{b}$ is fed into the trained DNN to output bootstrapping
samples $\hat{\bm{\theta}}_{b}$. The spread of $\hat{\bm{\theta}}_{b}$ and its
distance from $\hat{\bm{\theta}}_{0}$ are used to guide simulations for
training the DNN. Empirical results show a higher concentration of training
data and bootstrapping samples around the true parameters after a few
iterations, hence a better approximation of $\bm{\theta}$.

  2. (B)

A general unified database approach for time series: Computation efficiency is
achieved by training the DNN in advance and reusing it to estimate newly
collected data for free. The pre-training is based on an extensive database
comprising simulated time series data $\mathbf{x}$ and corresponding
parameters $\bm{\theta}$. Next, parameter estimates from new data of different
lengths are obtained by replicating the observations to achieve the size of
$\mathbf{x}$.

The details on the building blocks of the frameworks in (A) and (B), along
with the necessary adjustments for different applications, are given in
Sections 3 and 4.

##  3 A fully automatic iterative approach

###  3.1 General framework

We now use the notation in Section 2 to describe an algorithm that
sequentially samples training data until it reaches the correct parameter
regions from the data. Our algorithm is initialized with simulated pairs
($\bm{\theta}_{n},\mathbf{x}_{n})^{N}_{n=1}$, where the elements in
$\bm{\theta}_{n}\in\mathbb{R}^{P}$ are draw from $P$ independent Uniform
distributions, each bounded below by $(a_{1,p})^{P}_{p=1}$ and above by
$(a_{2,p})^{P}_{p=1}$, while $\mathbf{x}_{n}\in\mathbb{R}^{J}$ is data
simulated with $\bm{\theta}_{n},n=1,\ldots,N$. Intervals
$(a_{1,p},a_{2,p})^{P}_{p=1}$ may or may not contain the true parameter
$(\theta_{0,p})^{P}_{p=1}$. Whereas good initial guesses of
$(a_{1,p})^{P}_{p=1}$ and $(a_{2,p})^{P}_{p=1}$ are not essential here, most
models allow for some data-driven estimates, e.g., based on simplified
Gaussian assumptions. Next, a DNN is trained with
($\bm{\theta}_{n},\mathbf{x}_{n})^{N}_{n=1}$, and then used to retrieve
estimates $\hat{\bm{\theta}}_{0}\in\mathbb{R}^{P}$ when fed with observations
$\mathbf{x}_{0}\in\mathbb{R}^{J}$.

The main idea is to then dynamically update the training data
($\bm{\theta}_{n},\mathbf{x}_{n})^{N}_{n=1}$ by changing the values of
$(a_{1,p})^{P}_{p=1}$ and $(a_{2,p})^{P}_{p=1}$ based on information on
whether $\hat{\theta}_{0,p}$ is underestimating or overestimating
${\theta}_{0,p}$. For instance, if $\hat{\theta}_{0,p}$ is close to the upper
boundary of the training data for a specific $p$, then new training samples
should be expanded to contain data outside of that boundary. Information on
the accuracy of $(\hat{\theta}_{0,p})^{P}_{p=1}$ is obtained by sampling new
data $\mathbf{x}_{b}=(\mathbf{x}_{b}^{1},\ldots,\mathbf{x}_{b}^{B})^{\top}$
with $\mathbf{x}_{b}\in\mathbb{R}^{J\times B}$ using $\hat{\bm{\theta}}_{0}$,
and feeding these data into the initially trained DNN producing a bootstrapped
sample
$\hat{\bm{\theta}}_{b}=(\hat{\bm{\theta}}_{b}^{1},\ldots,\hat{\bm{\theta}}_{b}^{B})^{\top},\hat{\bm{\theta}}_{b}\in\mathbb{R}^{P\times
B}$. For each $p$, we then update $a_{1,p}$ and $a_{2,p}$ to values in a
neighborhood of $\hat{\theta}_{0,p}$, where the neighborhood region is defined
by the size of the bias between $\hat{\theta}_{0,p}$ and $\tilde{\theta}_{p}$,
where $\tilde{\theta}_{p}$ is the median of
$\theta^{1}_{p},\ldots,\theta^{B}_{p}$, and the neighborhood width depends on
the quantiles of the bias between the fitted value and the bootstrapped
sample:
$\mathcal{Q}^{\alpha}_{p}(\hat{\theta}_{0,p}-\theta^{1}_{p},\ldots,\hat{\theta}_{0,p}-\theta^{B}_{p})$,
where $\alpha$ is a quantile. The algorithm stops when the bias between
$\hat{\theta}_{0,p}$ and $\tilde{\theta}_{p}$ is sufficiently small compared
to the standard deviation of $(\theta^{1}_{p},\ldots,\theta^{B}_{p})^{\top}$,
which we denote by $S_{p}$, for all $p$. In more detail, the algorithm is as
follows.

Algorithm 1 Iterative procedure

Need: Observations $\mathbf{x}_{0}\in\mathbb{R}^{J}$ from a distribution
$p,\bm{\theta}_{0}\in\mathbb{R}^{P}$ and a neural network
$\mathcal{F}_{\phi}(\cdot)$
Pick $\gamma\in(0,1)$, $a_{1,p}$ and $a_{2,p},p=1,\ldots,P,$

1:while $\mbox{bias}(\hat{\theta}_{0,p},\tilde{\theta}_{b,p})>\gamma\times
S_{p}$, for all $p$ do

2: Sample $\theta_{n,p}\sim\mbox{Unif}(a_{1,p},a_{2,p}),n=1,\ldots,N$

3: Simulate $\mathbf{x}_{n}\sim p(;\bm{\theta}_{n}),n=1,\ldots,N$

4: Train $\mathcal{F}_{\phi}(\mathbf{x})$ and obtain $\hat{\bm{\theta}}_{0}$
from $\mathcal{F}_{\hat{\phi}}(\mathbf{x}_{0})$

5: Simulate $\mathbf{x}_{b}\sim p(;\hat{\bm{\theta}}_{0})$ and obtain
$\hat{\bm{\theta}}_{b}$ from
$\mathcal{F}_{\hat{\phi}}(\mathbf{x}_{b}),b=1,\ldots,B$

6:
$a_{1,p}=\hat{\theta}_{0,p}+\mbox{bias}(\hat{\theta}_{0,p},\tilde{\theta}_{p})-\mathcal{Q}^{0.05}_{p}(\hat{\theta}_{0,p}-\theta^{1}_{p},\ldots,\hat{\theta}_{0,p}-\theta^{B}_{p})$

7:
$a_{2,p}=\hat{\theta}_{0,p}+\mbox{bias}(\hat{\theta}_{0,p},\tilde{\theta}_{p})+\mathcal{Q}^{0.975}_{p}(\hat{\theta}_{0,p}-\theta^{1}_{p},\ldots,\hat{\theta}_{0,p}-\theta^{B}_{p})$

8:end while

Small values of $\gamma$ in line 1 of Algorithm 1 will make the algorithm run
longer, since it requires $\hat{\theta}_{0,p}$ to be closer to
$\tilde{\theta}_{b,p}$ relative to the spread of
$\theta^{1}_{p},\ldots,\theta^{B}_{p}$. At each iteration, line 5
automatically provides uncertainty quantification of $\hat{\bm{\theta}}_{0}$
through $\hat{\bm{\theta}}_{b}$. This step works as a modified and more
efficient parametric bootstrap method since it uses the previously trained
DNN, and no model fitting is required to produce $\hat{\bm{\theta}}_{b}$. One
can use these samples to compute quantities of interest, such as confidence
intervals and coverage, and check the overall appropriateness of the method.
Here, we use them to quantify the accuracy of the current iteration and to
update the training data for the next round (see lines 6 and 7).

Parameter values $\bm{\theta}$ are usually in the transformed scale, and we
continuously sample training data such that the values in this transformed
scale are uniformly distributed (see line 2 of Algorithm 1) rather than
applying transformations after the training data have been generated to train
the DNN. The former would produce regions of scarcity in $\Theta$, and results
for testing data within the underrepresented values would not be optimal.
Indeed, the optimization inside the DNN will perform best if the training data
have no significant gaps between values, a problem also called imbalanced data
in classification problems (Murphey et al.,, 2004).

In Section 3.2, we estimate the parameters of an i.i.d Gaussian model as a
proof-of-concept, whereas, in Section 3.3, we consider a spatial-extremes
setting and estimate the parameters of the Brown-Resnick max-stable process
with an intractable likelihood. This procedure supports a wide range of
likelihoods with fixed and random effects, and distributions other than the
Uniform could also have been used in line 2 of Algorithm 1.

###  3.2 I.i.d. data

Although classical inference for the models considered in this section is
straightforward, they allow us to compare our estimates’ accuracy and
uncertainty with MLEs. For applications where our method is of practical
interest, see Section 3.3. In what follows, we look at three problems of
increasing complexity from parameters of Gaussian distributions: the logarithm
variance (single parameter), the mean and the logarithm variance (two
orthogonal parameters), and the first moment and logarithm of the second
moment (two highly dependent parameters). With these examples, we aim to
empirically illustrate that our framework: (1) approaches the MLE even when
the initial training data is relatively far from the actual value, and (2)
reaches the truth quicker when using meaningful parametrizations and
orthogonal parameters.

Consider i.i.d. observations $\mathbf{x}_{0}\in\mathbb{R}^{J}$ from a Gaussian
distribution $\mathcal{N}(\mu,\sigma_{0}^{2})$. We find that $\gamma=0.3$ in
line 1 of Algorithm 1 is enough to provide good estimation accuracy without
overly increasing computational cost. Algorithm 2 shows the steps of our
procedure for the i.i.d. case (see Algorithm 1 for the general case), whereas
some practical aspects are discussed in what follows.

Algorithm 2 Iterative procedure for i.i.d. data

Need: Observations $\mathbf{x}_{0}\in\mathbb{R}^{J}$ from a distribution
$p,\bm{\theta}_{0}\in\mathbb{R}^{P}$ and a neural network
$\mathcal{F}_{\phi}(\cdot)$
Pick $a_{1,p}$ and $a_{2,p},p=1,2$

1:while $\mbox{bias}(\hat{\theta}_{0,p},\tilde{\theta}_{b,p})>0.3\times
S_{p}$, for all $p$ do

2: Sample $\theta_{n,p}\sim\mbox{Unif}(a_{1,p},a_{2,p}),n=1,\ldots,N$

3: Simulate $\mathbf{x}_{n}^{*}\sim p(;\bm{\theta}_{n}),n=1,\ldots,N$

4: Train $\mathcal{F}_{\phi}(\mathbf{x})$ and obtain $\hat{\bm{\theta}}_{0}$
from $\mathcal{F}_{\hat{\phi}}(\mathbf{x}_{0})$

5: Simulate $\mathbf{x}_{b}\sim p(;\hat{\bm{\theta}}_{0})$ and obtain
$\hat{\bm{\theta}}_{b}$ from
$\mathcal{F}_{\hat{\phi}}(\mathbf{x}_{b}),b=1,\ldots,B$

6:
$a_{1,p}=\hat{\theta}_{0,p}+\mbox{bias}(\hat{\theta}_{0,p},\tilde{\theta}_{p})-\mathcal{Q}^{0.05}_{p}(\hat{\theta}_{0,p}-\theta^{1}_{p},\ldots,\hat{\theta}_{0,p}-\theta^{B}_{p})$

7:
$a_{2,p}=\hat{\theta}_{0,p}+\mbox{bias}(\hat{\theta}_{0,p},\tilde{\theta}_{p})+\mathcal{Q}^{0.975}_{p}(\hat{\theta}_{0,p}-\theta^{1}_{p},\ldots,\hat{\theta}_{0,p}-\theta^{B}_{p})$

8: Increase $N$ by $5\%$

9:end while

#### Multi-layer perceptron (MLP)

A perceptron is a single neuron model, and MLPs are the classical type of
neural network comprised of one or more layers of several neurons. It takes 1D
vectors as the input and learns nonlinear relationships between inputs and
outputs, making it a suitable choice for our i.i.d. regression problem. Due to
the simplicity of this toy example, we find that a small MLP with a single
hidden layer and 50 hidden units is enough to near the mapping between data
and parameters. We consider an MLP for $\mathcal{F}_{\phi}$ taking output in
$\mathbb{R}^{P}$ and input in $\mathbb{R}^{J}$.

#### Progressively increasing training accuracy

Our algorithm uses fewer training samples when estimation uncertainty is
larger, at the beginning of the algorithm, and more samples towards the end,
when more precision is required. When the estimates are close to stabilizing,
and the uncertainty has decreased, the number of samples in the training data
is set to increase by $5\%$ (see line 8 of Algorithm 2).

#### Results for a single parameter

Figure 1 displays the results for estimating
$\theta_{0}\equiv\mbox{log}(\sigma_{0}^{2})=1$ when the mean is known using
Algorithm 2. We set $N=10000$ and uniformly generate training output samples
$\\{\theta_{n}\\}_{n=1}^{N}\sim{\mbox{Unif}}(-2,1)$ and corresponding inputs
$\\{\mathbf{x}_{n}\\}_{n=1}^{N}\sim N\\{1,\mbox{exp}(\theta_{n})\\}$. The grey
boxes in this figure are the training data at each iteration, whereas fitted
values are represented by the red line, with the blue boxes showing
bootstrapped estimates for $B=10000$. After five iterations, the algorithm
approaches the MLE (green dashed line) with low uncertainty (see narrow blue
boxes). To quantify the appropriateness of our method, we compare $95\%$
central intervals from the bootstrapping estimates with the same interval from
the empirical variance in the data. The $95\%$ interval provides adequate
uncertainty of the MLP estimates, with bootstrap interval on $(0.84,2.91)$
compared to $(0.83,2.95)$ for the intended coverage probability of the MLE.

Figure 1: Boxplots of training data (grey boxes), fitted values (red line),
and bootstrapped samples (blue boxes) as iterations progress for estimating
$\mbox{log}(\sigma_{0}^{2})=1$ from a zero-mean Gaussian distributed sample of
size $J=20$ using Algorithm 2. The horizontal dashed line corresponds to the
MLE. Training datasets were initially simulated in the uniform interval, with
$N=10000$ samples.

#### Results for two independent/dependent parameters

We now increase the problem’s complexity and evaluate the performance of
Algorithm 2 when two parameters are estimated jointly. We use the same test
data as in the single parameter estimation case (see Figure 1), that is, data
from a Gaussian distribution with $\mu_{0}=1$ and
$\mbox{log}(\sigma_{0}^{2})=1$. Specifically, we look at two cases: 1.
estimating the mean $\mu_{0}$ and log-variance $\mbox{log}(\sigma_{0}^{2})$
and 2. estimating the first moment $m_{1}=\mu_{0}$ and the logarithm of the
second moment $m_{2}=\mbox{log}(\mu_{0}^{2}+\sigma_{0}^{2})$. Whereas in case
1, the parameters are independent, the MLP has to learn the relation between
data and highly dependent parameters in case 2. In both cases, initial
training data does not contain the actual parameters, such that: 1.
$\mu_{n}\sim{\mbox{Unif}}(-0.5,0.5)$, and
$\mbox{log}(\sigma_{n}^{2})\sim{\mbox{Unif}}(-2,1)$, 2.
$m_{1,n}\sim{\mbox{Unif}}(-0.5,0.5)$ and
$m_{2,n}\sim{\mbox{Unif}}\\{-0.5^{2}+\mbox{exp}(-2),0.5^{2}+\mbox{exp}(1)\\},n=1,\ldots,N$.
Similarly to the single parameter estimation, we fix $N=10000$ and $J=20$.
Figure 2 displays the estimates for Case 1. (top row) and 2. (bottom row)
after running Algorithm 2 with the logarithm transformations (left column) and
without (right column). Estimates are more accurate when the parameters are
transformed, and for Case 2, the MLP underestimates the raw moments.
Especially for the second moment, the algorithm without transformation narrows
the estimates close to the median of the initial training data. In contrast,
the reparametrization is able to recover highly dependent parameters
accurately. Indeed, in both cases, the MLP can detect accurate parameter
estimates already in the first iteration. The subsequent iterations refine the
estimates and concentrate the training and the bootstrapping samples around
the MLEs. This experiment reiterates the benefit of using unbounded and
orthogonal parameters for training the MLPs.

(a) Mean and transformed variance

(b) Mean and variance

(c) First and second moment transformed

(d) First and second moments

Figure 2: Boxplots of the training data (grey) and bootstrap uncertainty
(blue) at different iterations of Algorithm 2. Points in the red line are the
fitted values from the MLP and the green dashed lines are the MLEs. Training
output data were initially simulated using $N=10000$ training samples each of
length $J=20$ as: (a) $\mu_{n}$ and log($\sigma_{n}^{2}$), (b) $\mu_{n}$ and
$\sigma_{n}^{2}$, (c) $\mu_{n}^{2}$ and log($\mu_{n}^{2}+\sigma_{n}^{2}$), and
(d) $\mu_{n}^{2}$ and $\mu_{n}^{2}+\sigma_{n}^{2}$.

###  3.3 Spatial extremes

We now move to a more complex model for spatial extremes, which is well-known
to have a likelihood function that is effectively impossible to compute. Max-
stable distributions are the only possible non-degenerate limits of
renormalized pointwise maxima of i.i.d random fields and, therefore, the most
commonly used for studying multivariate extreme events Davison et al., (2012).
We consider the following definition of a max-stable process

$X(\mathbf{s})=\max\limits_{i\geq
1}\xi_{i}W_{i}(\mathbf{s}),\mathbf{s}\in\mathcal{S},$ (2)

where $\\{\xi_{i}\\}_{i\geq 1}$ are points of a Poisson process on
$(0,\infty)$ with intensity $d\Lambda(\xi)=\xi^{-2}d\xi$. We consider the
Brown-Resnick model (Kabluchko et al.,, 2009), which arises when
$W_{i}(\mathbf{s})=\mbox{exp}\\{\epsilon_{i}(\mathbf{s})-\gamma_{i}(\mathbf{s})\\}$.
Each $\\{W_{i}\\}_{i\geq 1}$ is a nonnegative stochastic process with unit
mean, whereas $\epsilon_{i}(\mathbf{s})$ are copies of a zero-mean Gaussian
process with semivariogram
$\gamma(\mathbf{h})=(\lVert\mathbf{h}\rVert/\lambda)^{\nu}$, spatial
separation distance $\mathbf{h}$, range $\lambda>0$, smoothness $\nu\in(0,2]$
and such that
$\sigma^{2}(\mathbf{h})=\mbox{var}\\{\epsilon(\mathbf{h})\\}=2\gamma(\mathbf{h})$.

The cumulative distribution of $X(\mathbf{s})$ is

$p(X(\mathbf{s}_{1})\leq x_{1},\ldots,X(\mathbf{s}_{D})\leq
x_{D})=\mbox{exp}\\{-V(x_{1},\ldots,x_{D})\\},$

where
$V(x_{1},\ldots,x_{D})=\mbox{E}[\mbox{max}\\{W(\mathbf{s}_{1})/x_{1},\ldots,W(\mathbf{s}_{D})/x_{D}\\}]$
satisfies homogeneity and marginal constraints. The full likelihood is written
as

$f(x_{1},\ldots,x_{D})=\mbox{exp}\\{-V(x_{1},\ldots,x_{D})\\}\sum_{\pi\in\Gamma}\prod_{r=1}^{R}\\{-V_{\pi_{r}}(x_{1},\ldots,x_{D})\\},$

where $\Gamma$ is a collection of all partitions
$\pi=\\{\pi_{1},\ldots,\pi_{R}\\}$ of $\\{1,\ldots,D\\}$ and $V_{\pi_{r}}$
denotes the partial derivative of $V$ with respect to the variables indexed by
$\pi_{r}$. The full likelihood is intractable even for moderate $D$ since the
number of terms grows equals the Bell number, which is more than
exponentially. The standard workaround for this issue is to consider only
pairs of possibly weighted observations in the likelihood (Padoan et al.,,
2010; Davis et al.,, 2013; Shang et al.,, 2015):

$l(\bm{\theta})=\sum_{(j_{1},j_{2})\in\mathcal{P}}\alpha_{j_{1},j_{2}}\Big{[}\mbox{log}\\{V_{1}(x_{j_{1}},x_{j_{2}})V_{2}(x_{j_{1}},x_{j_{2}})-V_{12}(x_{j_{1}},x_{j_{2}})\\}-V_{1}(x_{j_{1}},x_{j_{2}})\Big{]},$

where $x_{j}$ is the block maximum at location $j$,
$\Gamma=\\{(j_{1},j_{2}):1\leq j_{1}<j_{2}\leq D\\}$,
$\bm{\theta}\in\bm{\Theta}\subset\mathbb{R}^{P}$ is the vector of unknown
parameters and $\alpha_{j_{1},j_{2}}\geq 0$ is the weight of
$\\{j_{1},j_{2}\\}$.

We compare the estimators from our fully automatic iterative approach with
pairwise likelihood estimation on $I=100$ independent simulated datasets of a
Brown-Resnick model with $\lambda_{0}=6.2$ and $\nu_{0}=1$ on a spatial domain
$\mathcal{S}$ of size $[0;30]\times[0;30]$ with unit-square grid cells. The
steps used for estimating parameters of Brown-Resnick processes are shown in
Algorithm 3. To improve accuracy and efficiency, the pairwise likelihood if
fit only with pairs with at most 5-units apart and using the R-function
fitmaxstab from the SpatialExtremes R-package (Ribatet,, 2013).

Algorithm 3 Iterative procedure for Brown-Resnick processes

Need: Observations $\mathbf{x}_{0}\in\mathbb{R}^{J}$ from a distribution
$p,\bm{\theta}_{0}\in\mathbb{R}^{P}$ and a neural network
$\mathcal{F}_{\phi}(\cdot)$
Pick $a_{1,p}$ and $a_{2,p},p=1,2$, and set $\mathcal{D}=\\{\\}$

1:while $\mbox{bias}(\hat{\theta}_{0,p},\tilde{\theta}_{b,p})>0.3\times
S_{p}$, for all $p$ do

2: Sample $\theta_{n,p}\sim\mbox{Unif}(a_{1,p},a_{2,p}),n=1,\ldots,N$

3: Simulate $\mathbf{x}_{n}\sim p(;\bm{\theta}_{n}),n=1,\ldots,N$

4: Set
$\mathcal{D}_{\mbox{\scriptsize{train}}}=({\bm{\theta}}_{n},\mathbf{x}_{n})^{N}_{n=1}$

5: Train $\mathcal{F}_{\phi}(\mathbf{x})$ on
$\mathcal{D}_{\mbox{\scriptsize{train}}}\bigcup\mathcal{D}$ and obtain
$\hat{\bm{\theta}}_{0}$ from $\mathcal{F}_{\hat{\phi}}(\mathbf{x}_{0})$

6: Simulate $\mathbf{x}_{b}\sim p(;\hat{\bm{\theta}}_{0})$ and obtain
$\hat{\bm{\theta}}_{b}$ from
$\mathcal{F}_{\hat{\phi}}(\mathbf{x}_{b}),b=1,\ldots,B$

7:
$a_{1,p}=\hat{\theta}_{0,p}+\mbox{bias}(\hat{\theta}_{0,p},\tilde{\theta}_{p})-\mathcal{Q}^{0.05}_{p}(\hat{\theta}_{0,p}-\theta^{1}_{p},\ldots,\hat{\theta}_{0,p}-\theta^{B}_{p})$

8:
$a_{2,p}=\hat{\theta}_{0,p}+\mbox{bias}(\hat{\theta}_{0,p},\tilde{\theta}_{p})+\mathcal{Q}^{0.975}_{p}(\hat{\theta}_{0,p}-\theta^{1}_{p},\ldots,\hat{\theta}_{0,p}-\theta^{B}_{p})$

9: Randomly select a subset of $\mathcal{D}_{\mbox{\scriptsize{train}}}$ such
that $\mathcal{D}_{\mbox{\scriptsize{train}}}\cap\mathcal{D}=\emptyset$ and
add those into $\mathcal{D}$

10:end while

#### Convolution neural network (CNN)

CNN uses convolutions, that is, the application of a filter to the input image
that results in what is called an activation. Repeated application of the same
filter to images results in a map of activations (feature map). This map
indicates the locations and strength of a detected feature in the input, such
as the edges of objects, and therefore, it is a common choice for regularly-
spaced gridded images. We use two 2D convolutions with 16 and 8 filters,
respectively, and rectified linear unit (ReLU) activation function and kernel
of size $3\times 3$ (Hastie et al.,, 2009). We add one dense layer at the end
of the network with four units that map the input image to an output vector of
size two. The CNN weights are initialized randomly and trained using the Adam
optimizer (Kingma and Ba,, 2014) with a learning rate of 0.01. The training is
performed for 30 epochs, and at each epoch, the CNN weights are updated
utilizing a batch size of 100 samples from the entire training dataset.

#### Initialization

We start the by simulating $N=6000$ pairs
$(\bm{\theta}_{n},\mathbf{x}_{n})^{N}_{n=1}$, where
$\bm{\theta}_{n}\equiv(\theta_{n,1},\theta_{n,2})^{\top}\equiv\\{\mbox{log}(\lambda_{n}),\mbox{logit}(\nu_{n})\\}^{\top}$
and $\mathbf{x}_{n}$ is simulated from (2). We initialize
$(\theta_{n,1})_{n=1}^{N}$ based on estimates of a Gaussian process with
powered exponential covariance function
$C(\mathbf{h})=\mbox{exp}(-(\lVert\mathbf{h}\rVert/\alpha)^{\eta}),\alpha>0,0<\eta\leq
2$, which closely matches the Brown-Resnick variogram. We sample
$\theta_{n,1}\sim\mbox{Unif}\\{\mbox{log}(\hat{\alpha_{n}})-c,\mbox{log}(\hat{\alpha_{n}})+c\\}$.
Such estimates are likely biased for Brown-Resnick, but they are quick to
compute and a better start than a random guess. A line search for $c$
indicates that $c=2$ provides good results but found that other values of $c$
gave similar accuracy. We uniformly draw $(\theta_{n,2})_{n=1}^{N}$ over
approximately the whole (bounded) domain:
$\theta_{n,2}\sim\mbox{Unif}\\{\mbox{logit}(0.1),\mbox{logit}(1.9)\\}$. We
initialize the pairwise likelihood with the powered exponential covariance
function estimates for a fair comparison with our approach.

If the parameter region for simulating data at each iteration is tightened too
quickly or shifted by too much, it may miss the true parameters. To alleviate
this issue, we randomly select $40\%$ of the samples from the previous
iteration that were not contained in the updated interval and add them to the
training data of the current step (see line 9 of Algorithm 3). Besides
increasing the size of the training without having to simulate new data, this
broadens the range of the training while still keeping most of the samples in
the updated region defined by the current step. Therefore, it prevents the
current iteration from being stuck in the wrong region while maintaining more
accurate training data in the most probable parameter region.

Figure 3 shows a scatterplot of 100 independent estimates of
$\theta_{1}\equiv\mbox{log}(\lambda_{0})$ versus
$\theta_{2}\equiv\mbox{log}(\nu_{0})$ from the last iteration of our approach
(green) and from the pairwise likelihood (red). The $\times$ symbol is the
truth. Whereas the proposed method produces robust results across the
different replicates, the pairwise likelihood tends to underestimate the
smoothness parameter, and the performance varies considerably across datasets.
In Figure 4, we access the accuracy of Algorithm 3 for estimating $\theta_{1}$
(left column) and $\theta_{2}$ (right column) with boxplots. The rows in this
figure illustrate the results for two different datasets: The first is
initialized with training data that do not contain $\theta_{1}$ (top), whereas
the variogram estimate for the second dataset is close to the center of the
training data. Indeed, whereas the space covering $\theta_{2}$ is bounded, and
simulating training data covering the entire region is straightforward, the
training data for $\theta_{2}$ is based on the variogram estimate and,
therefore, only sometimes contains the truth. The grey and blue boxplots at
each plot and iteration are the training output and bootstrapping samples,
respectively. As expected, the gray boxplots in all cases at iterations 2 and
3 contain several outliers, which correspond to the samples reused from the
previous step. Points in the red line are the fitted values from the CNN, and
the green dashed lines are the actual parameters used to simulate data.

Even when the true value is not included in the initial training data (see the
top of Figure 4), the CNN estimates $\theta_{1}$ well and produces reasonable
uncertainties. At the last iteration, the training data for both parameters
are narrow and around the actual value, and the bootstrapping samples
practically coincide with the training data. When initialized with training
data containing the truth, the CNN quickly approaches the truth with low
uncertainty for both parameters and remains stable until it reaches the
stopping criteria. A quantitative measure of the effect of the iterations in
Algorithm 3 is reported in Table 1, with bias, standard deviation, and root
mean square error (RMSE) for the first and last iterations. The metrics are
calculated from the bootstrapping samples among the 100 independent datasets.
Point estimates are taken as the median of the bootstrapping samples. Under
all three metrics, there is a considerable improvement from the first to the
last iteration of the algorithm, and estimation at the last iteration are
about $15\%$ more efficient for $\theta_{1}$ and $30\%$ more efficient for
$\theta_{2}$ (with the efficiency defined as the ratio of RMSEs).

Figure 3: Scatterplots of estimated parameters on the transformed scales. Each dot/cross shows 100 independent estimates from the Brown-Resnick model using the CNN (green) or PL (red). The $\times$ is the actual value. Training datasets were initially simulated using $N=6000$ training samples on a $[0,30]^{2}$ and based on estimates from fitting a powered exponential covariance function to the data. Figure 4: Boxplots of the training data (grey) and bootstrap uncertainty (blue) at different iterations of Algorithm 3 for $\theta_{1}$ (left) and $\theta_{2}$ (right). The top/bottom row shows an example of when the variogram estimates are outside/inside the initial training data. Points in the red line are the fitted values from the CNN, and the green dashed lines are the truth. | $\theta_{1}$ | $\theta_{2}$
---|---|---
Iteration | bias | sd | RMSE | bias | sd | RMSE
First | 0.223 | 0.295 | 0.563 | 0.005 | 0.253 | 0.451
Last | 0.056 | 0.129 | 0.306 | 0.004 | 0.120 | 0.313
Table 1: Standard deviation, bias and RMSE of the estimated parameters from
the first and last iteration of Algorithm 3

##  4 A general unified database approach for time series

###  4.1 General framework

Suppose we observe time series data
$\mathbf{x}_{0}=\\{x_{0}(1),\ldots,x_{0}(T)\\}^{\top}$ from a strictly
stationary process $\\{X_{0}(t):t\in\mathcal{T}\\}$ indexed on the temporal
domain $T\subset\mathbb{R}_{+}$. Let $p(·;\bm{\theta}_{0})$ be the probability
distribution of $X_{0}(t)$ depending on the parameter set
$\bm{\theta}_{0}\in\Theta\subset\mathbb{R}^{P}$. The stationarity assumption
is that the joint probability distribution of
$\\{X_{0}(t-l),\ldots,X_{0}(t)\\}^{\top}$ does not depend on
$X_{0}(t-l^{\prime})$ for any $l^{\prime}>l$. This Markov property is common
in time series analysis similarly with ergodicity, which provides the
justification for estimating $\bm{\theta}_{0}$ from a single sequence
$\mathbf{x}_{0}$.

Our main goal is to estimate $\bm{\theta}_{0}$ by training a DNN
$\mathcal{F}_{\phi}$ using parameter candidates as output and corresponding
simulated data as input (see Section 2). Here, we take advantage of the
stationarity property to generalize and improve the estimation workflow
described in Section 2. Our approach is best exemplified by a toy data
$\mathbf{x}_{0}=\\{x_{0}(1),\ldots,x_{0}(T)\\}^{\top}$ with $T=50$ from an AR
process of order 1 with coefficient $\rho_{1}=0.9$. Instead of simulating
training data of length $T=50$, we proceed by simulating time series of length
$T_{k}=250$ and construct data pairs
$(\bm{\theta}_{n},\mathbf{x}_{n})^{N}_{n=1}$, where
$\bm{\theta}_{n}\in\mathbb{R}^{P}$ and
$\mathbf{x}_{n}=\\{x_{n}(1),\ldots,x_{n}(250)\\}^{\top}$ and train a DNN.
Next, since $\mathbf{x}_{0}$ is shorter than the training data, we create a
new time series $\mathbf{x}^{*}_{0}$ by concatenating $\mathbf{x}_{0}$ to
achieve the length of the time series used during training. Figure 5 shows how
to construct $\mathbf{x}^{*}_{0}$ from $\mathbf{x}_{0}$ by replicating the
observations five times. The red dashed line are the joining points. The
resulting $\mathbf{x}^{*}_{0}$ is then fed into the trained DNN to retrieve
estimates $\hat{\bm{\theta}}_{0}$.

Figure 5: Data from an autoregressive process of order one (AR(1)) of length
$T=50$ with AR coefficient equal to 0.9, augmented five times to achieve the
training data length $T_{k}=250$.

Recurrent neural networks take information from prior inputs to predict the
current input and output and can handle sequence data and inputs of varying
lengths by memorizing historical information. Therefore, these networks are
unsuitable for our problem, as our output data does not have future sequence
values. Here, we are interested in mapping parameters to sequential data
rather than learning the next value in a sequence. As we will see in the
example in the next section, our technique has several advantages over
training the DNN using simulated data with the same length as the
observations: (i) allows estimation of time series of several sizes at almost
no computational cost, since the DNN does not have to be retrained for each
new dataset, (ii) holds without requiring any particular dependence structure
assumption as long as the data is stationary and Markov.

#### Connection to non-overlapping block bootstrap (NBB)

The intuition behind our approach resembles NBB approaches (Carlstein,, 1986).
Similarly, this technique splits the observations into non-overlapping blocks
and resamples the blocks with replacement, which are replicated to obtain a
bootstrapped series. However, unlike block bootstrap methods, where the
complex problem of choosing the block size has to be solved, by construction,
our block size is always fixed and equal to $T$.

#### Discontinuity at the joining points

Our procedure of laying sequences of length $T$ end-to-end will inevitably
produce $m-1$ discontinuity points where the joining occurs, similarly to what
happens in bootstrap for time series. However, as we will show empirically in
our examples, these discontinuities will have a negligible contribution to the
model parameters structure.

#### Varying observation lengths

As long as the database for training the neural network is extensive enough,
the proposed method is easily generalized for cases where $T$ is not a
multiple of $T_{k}$. If the observed time series is shorter than that used
during training, that is $T<T_{k}$, one can replicate the data into $m$
blocks, where $m=\lfloor\frac{T_{k}}{T}\rfloor$, and complete the remaining
values with a random block from $\mathbf{x}_{0}$ of size $T_{k}\mbox{mod}T$.
The idea is that if $T_{k}$ is large enough compared to $T$ and with
$mT\backsimeq T_{k}$, as $m\rightarrow\infty$ and $T_{k}\rightarrow\infty$,
the last components of $\mathbf{x}^{*}_{0}$ have little influence in the
dependence structure. On the other hand, if $T>T_{k}$, we take the first
$T_{k}$ observations and estimate parameters $\bm{\theta}^{(1)}_{0}$, the next
$T_{k}$ and estimate $\bm{\theta}^{(2)}_{0}$, and so on. If the last sequence
from this procedure is smaller than $T_{k}$, we use the method for dealing
with $T<T_{k}$ described previously. Next, we merge the results from these
several estimates (e.g., by using the sample average) to create a single
combined estimate $\hat{\bm{\theta}}_{0}$, which can then be used to simulate
bootstrap samples and construct confidence intervals.

#### Uncertainty quantification

To account for the underestimation in the uncertainty that results from
applying the estimator to the replicated time series, we scale the sampling
variance by the number of blocks used to replicate the original series. For
instance, consider the simulated toy example in Figure 5 of an AR(1) process
with $T=200$ and an AR coefficient $\rho=0.9$. The MLE is
$\hat{\rho}_{T}\approx 0.8320$, and its corresponding standard deviation is
$0.0784$. If the series is then replicated $5$ times such that the new length
is $T_{k}=250$, the new MLE is not affected ($\hat{\rho}_{T_{k}}\approx
0.8320$), but its standard deviation is now $0.0357$. Therefore, to recover
the correct uncertainty, we need to scale the new standard deviation by
$\sqrt{5}$, that is, $0.0357\times\sqrt{5}\approx 0.0798$, which is then,
approximately, the standard deviation of $\hat{\rho}_{T}$. This adjusts the
uncertainty for the fact that no additional information is added to the
replicated data.

#### 1D Convolutional Neural Networks

Since we are dealing with parameters from time series data, we train 1D CNNs,
which have proven successful in learning features from dependent observations
onto one-dimensional dependent sequences. As for 2D CNNs (see Section 3.3 for
an example), the input layer receives the (transformed) data, and the output
layer is an MLP with the number of neurons equal to the number of output
variables. Each neuron in a hidden layer first performs a sequence of
convolutions, the sum of which is passed through the activation function
followed by a sub-sampling operation. The early convolutional layers can be
seen as smoothing the input vector, where the filters are similar to
parameters of a weighted moving average but learned jointly with the
regression parameters from the MLP layer. In what follows, we use three 1D
convolutions with four filters each, the ReLU activation function, a kernel of
size three, and one final dense layer with four units. We set a learning rate
of 0.01 with 30 epochs and a batch size of 50 samples to update the weights.

###  4.2 Non-Gaussian stochastic volatility model

To show the usefulness of our approach, we focus on estimating parameters of
financial time series data that exhibit non-Gaussian time-varying volatility.
Volatility is highly right-skewed and bounded, making Gaussian distributions a
poor representation. A better description of volatility is achieved with
stochastic volatility models (SVOL), first introduced in Taylor, (1982) and
currently central to econometrics and finance investments theory and practice.
The idea of SVOL models is to parsimoniously fit the volatility process as a
latent structure using an unconditional approach that does not depend on
observations. We consider the following model structure

$\displaystyle x(t)$
$\displaystyle=\sqrt{\mbox{exp}\\{h(t)\\}}\epsilon_{t},\quad\sqrt{\frac{\nu-2}{\nu}}\epsilon(t)\sim
T_{\nu},\quad t=1,\ldots,T$ (3) $\displaystyle h(t)$ $\displaystyle=\rho
h(t-1)+\xi_{t},\quad\xi(t)\sim N(0,\sigma^{2}),$

where $\epsilon(t)$ and $\xi(t)$ are independent noises and $T$ is the number
of observations. The volatility variable $h(t)$ is latent with an AR(1)
structure, and only $x(t)$ is observed. When $\lvert{\rho}\rvert<1$, $x(t)$ is
strictly stationary with mean $\mu_{h}=1/(1-\rho)$ and variance
$\sigma^{2}_{h}=\sigma^{2}/(1-\rho^{2})$ (Fridman and Harris,, 1998).

Likelihood evaluation of continuous dynamic latent-variable models such as (3)
requires the integration of the latent process out of the joint density,
resulting in the following $T$-dimensional integral

$\displaystyle l({\bm{\theta}})$
$\displaystyle=\int_{\mathbf{h}}p(\mathbf{x}\mid\mathbf{h},{\bm{\theta}})p(\mathbf{h}\mid{\bm{\theta}})d\mathbf{h}$
(4)
$\displaystyle=\int_{\mathbf{h}}\prod^{T}_{t=1}p(\mathbf{x}\mid\mathbf{h})p\\{h(t)\mid
h(t-1)\\}d\mathbf{h},$

where $\mathbf{x}=\\{x(1),\ldots,x(T)\\}^{\top}$,
$\mathbf{h}=\\{h(1),\ldots,h(T)\\}^{\top}$, the conditional densities
$p(\mathbf{x}\mid\mathbf{h},{\bm{\theta}})$ and $p\\{h(t)\mid h(t-1)\\}$ have
the form in (3), and the initial volatility $h(0)$ is the stationary
volatility distribution $p(\mathbf{h})$. Since $h(t)$ is not independent from
the past, the integral in (4) cannot be factored into a product of $T$ one-
dimensional integrals and exact evaluation of the likelihood is possible only
in special cases like when both $p(\mathbf{x}\mid\mathbf{h},{\bm{\theta}})$
and $p\\{h(t)\mid h(t-1)\\}$ are Gaussian. Alternative approaches for
likelihood evaluation include computationally demanding Markov Chain Monte
Carlo (MCMC) (Andersen et al.,, 1999) and the more recent and faster
Integrated Nested Laplace Approximation (INLA) (Martino et al.,, 2011).

Next, we give practical details of our framework as well as the a comparison
of the results from our approach and the state-of-the-art INLA approach for
estimating parameters of the SVOL model.

#### Implementation

Consider observations $\mathbf{x}_{0}=\\{x_{0}(1),\ldots,x_{0}(T)\\}^{\top}$,
from the SVOL model with $\sigma_{0}=0.1$, $\rho_{0}=0.8$ and $\nu_{0}=6$. We
use scaled (variance one) versions of both $h(t)$ and $x(t)$ in (3), such that
only $\rho_{0}$ and $\nu_{0}$ need to be estimated. To show the effect of
estimating time series of different lengths with a single DNN fit, we display
the results for various time series lengths:
$T=(500,1000,2000,3000,4000,5000)$. Estimation goes as follows. The training
database contain $N=10000$ samples pairs
$(\bm{\theta}_{n},\mathbf{x}_{n})^{10000}_{n=1}$ of transformed parameters
$\bm{\theta}_{n}=(\theta_{1},\theta_{2})^{\top}=\\{f_{1}(\rho_{n}),f_{2}(\nu_{n})\\}^{\top}$,
with $f_{1}(x)=\mbox{log}\Big{(}\frac{1+x}{1-x}\Big{)}$ and
$f_{2}(x)=\mbox{log}(x-2)$ and corresponding data
$\mathbf{x}_{n}=\\{x_{n}(1),\ldots,x_{n}(5000)\\}^{\top}$ simulated from (3).
The transformed parameters are sampled uniformly in a neighborhood of the
actual parameter values $\theta_{0,1}=f_{1}(\rho_{0})$ and
$\theta_{0,2}=f_{2}(\nu_{0})$:

$\displaystyle f_{1}(\rho_{j})$
$\displaystyle\sim{\mbox{Unif}}\\{\theta_{0,1}-c^{a}_{1},\theta_{0,1}+c^{b}_{1}\\}$
(5) $\displaystyle f_{2}(\nu_{j})$
$\displaystyle\sim{\mbox{Unif}}(\theta_{0,2}-c^{a}_{2},\theta_{0,2}+c^{b}_{2}).$

New test data $\mathbf{x}^{*}_{0}$ is obtained by replicating $\mathbf{x}_{0}$
$m$ times as many time as needed to achieve size $5000$. We fix
$c^{a}_{1}=c^{a}_{2}=c^{b}_{1}=c^{b}_{2}=2$ in (LABEL:eq:ts_sim) to ensure a
large enough region around the true values, although other constants provided
similar results.

Figure 6 displays scatterplots of estimated $\theta_{0,1}$ versus estimated
$\theta_{0,2}$ from $I=30$ independent replicates of the SVOL model.
Scatterplots from top to bottom, left to right, shows testing sets of
different size: $T=(500,1000,2000,3000,4000,5000)$. Green dots are the 1D CNN
estimates, and red dots are the mean of the predicted posterior distribution
from fitting model (3) using INLA. The $\times$ symbol in each plot represents
the truth. As the testing data sizes increase, both methods concentrate the
estimates around the truth, and the 1D CNN estimates are less variable and
less biased for smaller values of $T$; after that, both methods seem to
perform similarly well. Whereas we set the INLA priors to default values,
changing them to penalize parameter values far from the mode could potentially
improve the results.

To quantify the uncertainty in the 1D CNN estimates, we create a bootstrapped
dataset by independently sampling $B=100000$ time series from the fitted model
and then feeding these samples into the trained 1D CNN. The uncertainty in the
estimation of $\theta_{0,1}$ (left) and $\theta_{0,2}$ (right) from both
methods is grasped in Figure 7. The bars represent central $95\%$ intervals,
taken from the posterior distribution given by INLA (red) and the 1D CNN
bootstrapped samples (blue) for one randomly chosen dataset among the $I=30$
replicates in Figure 6. The $x$-axis represents different test data sizes
$T=(500,1000,2000,3000,4000,5000)$, and the gray horizontal dashed line is the
truth. As expected, the uncertainties decrease with sample sizes for both
methods. The intervals from both methods are close for most data sizes,
although the INLA distributions are more concentrated for $\theta_{0,1}$ and
$T=(4000,5000)$. Overall, these results show that the 1D CNN is robust to
estimating parameters of different data lengths and in agreement with the INLA
estimator.

Figure 6: Scatterplots of estimated parameters on the transformed scales. Each
plot shows 30 independent estimates from the SVOL model from the 1D CNN
(green) and INLA (red). Small testing data to large are displayed from top to
bottom, left to right: $T=\\{500,1000,2000,3000,4000,5000\\}$. The $\times$’s
are the true values.

Figure 7: Example of $95\%$ estimated central interval for $\theta_{0,1}$
(left) and $\theta_{0,2}$ (right) using the estimated posterior distribution
from INLA (red) and 10000 bootstrapped samples from the 1D CNN (blue). The
$x$-axis are the different test data sizes $T=(500,1000,2000,3000,4000,5000)$
and the horizontal dashed line is the truth.

##  5 Conclusion

We proposed approaches that train DNNs to estimate parameters of intractable
models and quantify their uncertainty. Unlike previously proposed approaches
using DNNs, which are tailored to a specific application and can lead to poor
parameter estimates for relying on computationally expensive initial guesses
to construct training data, our methods (A) leverage an iterative learning
framework coupled with a modified parametric bootstrap step to guide
simulations in the direction of the parameter region of the actual data in
multiple rounds (B) use an extensive pre-trained database to accurately
estimate parameters of time series data of multiple lengths at no
computational cost, rather than simulating data for every new dataset.

Our estimators yield accurate parameter estimates with much less computation
time than classical methods, even when accounting for the time required to
generate training samples. The experiments discussed in this work involve
models with a maximum of two parameters for clarity in presentation. We
anticipate a more pronounced influence of our sequential approach in higher-
dimensional scenarios for alleviating the challenges associated with the
“curse of dimensionality” as the volume of the parameter space increases
exponentially with the number of parameters.

While DNNs for parameter estimation are gaining popularity, we still need to
learn more about black-box algorithms applied to previously intractable
statistical problems and how to design task-specific estimators more
generally. There are several further opportunities for exploring DNNs for
parameter estimation using newly designed optimization tools from the machine
learning community. This work is another step towards this direction, where
ultimately, inference is performed within a general and flexible simulation-
based black-box pipeline.

## References

  * Andersen et al., (1999) Andersen, T. G., Chung, H.-J., and Sørensen, B. E. (1999).  Efficient method of moments estimation of a stochastic volatility model: A monte carlo study.  Journal of econometrics, 91(1):61–87.
  * Beaumont et al., (2009) Beaumont, M. A., Cornuet, J.-M., Marin, J.-M., and Robert, C. P. (2009).  Adaptive approximate bayesian computation.  Biometrika, 96(4):983–990.
  * Bonassi and West, (2015) Bonassi, F. V. and West, M. (2015).  Sequential monte carlo with adaptive weights for approximate bayesian computation.
  * Carlstein, (1986) Carlstein, E. (1986).  The use of subseries values for estimating the variance of a general statistic from a stationary sequence.  The annals of statistics, pages 1171–1179.
  * Cranmer et al., (2020) Cranmer, K., Brehmer, J., and Louppe, G. (2020).  The frontier of simulation-based inference.  Proceedings of the National Academy of Sciences, 117(48):30055–30062.
  * Davis et al., (2013) Davis, R. A., Klüppelberg, C., and Steinkohl, C. (2013).  Statistical inference for max-stable processes in space and time.  Journal of the Royal Statistical Society: SERIES B: Statistical Methodology, pages 791–819.
  * Davison et al., (2012) Davison, A. C., Padoan, S. A., and Ribatet, M. (2012).  Statistical modeling of spatial extremes.  Statistical science, 27(2):161–186.
  * Drovandi and Frazier, (2022) Drovandi, C. and Frazier, D. T. (2022).  A comparison of likelihood-free methods with and without summary statistics.  Statistics and Computing, 32(3):42.
  * Fearnhead and Prangle, (2012) Fearnhead, P. and Prangle, D. (2012).  Constructing summary statistics for approximate bayesian computation: semi-automatic approximate bayesian computation.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 74(3):419–474.
  * Frazier et al., (2018) Frazier, D. T., Martin, G. M., Robert, C. P., and Rousseau, J. (2018).  Asymptotic properties of approximate bayesian computation.  Biometrika, 105(3):593–607.
  * Fridman and Harris, (1998) Fridman, M. and Harris, L. (1998).  A maximum likelihood approach for non-gaussian stochastic volatility models.  Journal of Business & Economic Statistics, 16(3):284–291.
  * Friedman et al., (2001) Friedman, J., Hastie, T., Tibshirani, R., et al. (2001).  The Elements of Statistical Learning, volume 1.  Springer Series in Statistics.
  * Gerber and Nychka, (2020) Gerber, F. and Nychka, D. W. (2020).  Fast covariance parameter estimation of spatial Gaussian process models using neural networks.  Stat, page e382.
  * Gourieroux et al., (1993) Gourieroux, C., Monfort, A., and Renault, E. (1993).  Indirect inference.  Journal of applied econometrics, 8(S1):S85–S118.
  * Grelaud et al., (2009) Grelaud, A., Marin, J.-M., Robert, C. P., Rodolphe, F., and Taly, J.-F. (2009).  Abc likelihood-free methods for model choice in gibbs random fields.  Bayesian Analysis, 4(2):317–335.
  * Gutmann and Corander, (2016) Gutmann, M. U. and Corander, J. (2016).  Bayesian optimization for likelihood-free inference of simulator-based statistical models.  Journal of Machine Learning Research.
  * Hartig et al., (2011) Hartig, F., Calabrese, J. M., Reineking, B., Wiegand, T., and Huth, A. (2011).  Statistical inference for stochastic simulation models–theory and application.  Ecology letters, 14(8):816–827.
  * Hastie et al., (2009) Hastie, T., Tibshirani, R., and Friedman, J. (2009).  The Elements of Statistical Learning: Data Mining, Inference, and Prediction.  Springer Science & Business Media.
  * Järvenpää et al., (2019) Järvenpää, M., Gutmann, M. U., Pleska, A., Vehtari, A., and Marttinen, P. (2019).  Efficient acquisition rules for model-based approximate bayesian computation.
  * Jiang et al., (2017) Jiang, B., Wu, T.-y., Zheng, C., and Wong, W. H. (2017).  Learning summary statistic for approximate Bayesian computation via deep neural network.  Statistica Sinica, pages 1595–1618.
  * Kabluchko et al., (2009) Kabluchko, Z., Schlather, M., De Haan, L., et al. (2009).  Stationary max-stable fields associated to negative definite functions.  The Annals of Probability, 37(5):2042–2065.
  * Kingma and Ba, (2014) Kingma, D. P. and Ba, J. (2014).  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980.
  * Lenzi et al., (2021) Lenzi, A., Bessac, J., Rudi, J., and Stein, M. L. (2021).  Neural networks for parameter estimation in intractable models.  arXiv preprint arXiv:2107.14346.
  * Lueckmann et al., (2019) Lueckmann, J.-M., Bassetto, G., Karaletsos, T., and Macke, J. H. (2019).  Likelihood-free inference with emulator networks.  In Symposium on Advances in Approximate Bayesian Inference, pages 32–53. PMLR.
  * Lueckmann et al., (2017) Lueckmann, J.-M., Goncalves, P. J., Bassetto, G., Öcal, K., Nonnenmacher, M., and Macke, J. H. (2017).  Flexible statistical inference for mechanistic models of neural dynamics.  Advances in neural information processing systems, 30.
  * Martino et al., (2011) Martino, S., Aas, K., Lindqvist, O., Neef, L. R., and Rue, H. (2011).  Estimating stochastic volatility models using integrated nested laplace approximations.  The European Journal of Finance, 17(7):487–503.
  * Murphey et al., (2004) Murphey, Y. L., Guo, H., and Feldkamp, L. A. (2004).  Neural learning from unbalanced data.  Applied Intelligence, 21(2):117–128.
  * Nickl and Pötscher, (2010) Nickl, R. and Pötscher, B. M. (2010).  Efficient simulation-based minimum distance estimation and indirect inference.  Mathematical methods of statistics, 19(4):327–364.
  * Padoan et al., (2010) Padoan, S. A., Ribatet, M., and Sisson, S. A. (2010).  Likelihood-based inference for max-stable processes.  Journal of the American Statistical Association, 105(489):263–277.
  * Papamakarios and Murray, (2016) Papamakarios, G. and Murray, I. (2016).  Fast $\varepsilon$-free inference of simulation models with bayesian conditional density estimation.  In Advances in neural information processing systems, pages 1028–1036.
  * Papamakarios et al., (2019) Papamakarios, G., Sterratt, D., and Murray, I. (2019).  Sequential neural likelihood: Fast likelihood-free inference with autoregressive flows.  In The 22nd International Conference on Artificial Intelligence and Statistics, pages 837–848. PMLR.
  * Ribatet, (2013) Ribatet, M. (2013).  Spatial extremes: Max-stable processes at work.  Journal de la Société Française de Statistique, 154(2):156–177.
  * Richards et al., (2023) Richards, J., Sainsbury-Dale, M., Huser, R., and Zammit-Mangion, A. (2023).  Likelihood-free neural bayes estimators for censored peaks-over-threshold models.  arXiv preprint arXiv:2306.15642 (2023).
  * Rue et al., (2009) Rue, H., Martino, S., and Chopin, N. (2009).  Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 71(2):319–392.
  * Sainsbury-Dale et al., (2023) Sainsbury-Dale, M., Richards, J., Zammit-Mangion, A., and Huser, R. (2023).  Neural bayes estimators for irregular spatial data using graph neural networks.  arXiv preprint arXiv:2310.02600.
  * Sainsbury-Dale et al., (2022) Sainsbury-Dale, M., Zammit-Mangion, A., and Huser, R. (2022).  Fast optimal estimation with intractable models using permutation-invariant neural networks.  arXiv preprint arXiv:2208.12942.
  * Shang et al., (2015) Shang, H., Yan, J., and Zhang, X. (2015).  A two-step approach to model precipitation extremes in California based on max-stable and marginal point processes.  The Annals of Applied Statistics, pages 452–473.
  * Sisson et al., (2018) Sisson, S. A., Fan, Y., and Beaumont, M. (2018).  Handbook of approximate Bayesian computation.  CRC Press.
  * Sisson et al., (2007) Sisson, S. A., Fan, Y., and Tanaka, M. M. (2007).  Sequential monte carlo without likelihoods.  Proceedings of the National Academy of Sciences, 104(6):1760–1765.
  * Taylor, (1982) Taylor, S. J. (1982).  Financial returns modelled by the product of two stochastic processes-a study of the daily sugar prices 1961-75.  Time series analysis: theory and practice, 1:203–226.
  * Wood, (2010) Wood, S. N. (2010).  Statistical inference for noisy nonlinear ecological dynamic systems.  Nature, 466(7310):1102–1104.