Search

The Active Subspace (AS) method is a widely used technique for identifying the most influential directions in high-dimensional input spaces that affect the output of a computational model. The standard AS algorithm requires a sufficient number of gradient evaluations (samples) of the input output map to achieve quasi-optimal reconstruction of the active subspace, which can lead to a significant computational cost if the samples include numerical discretization errors which have to be kept sufficiently small. To address this issue, we propose a multilevel version of the Active Subspace method (MLAS) that utilizes samples computed with different accuracies and yields different active subspaces across accuracy levels, which can match the accuracy of single-level AS with reduced computational cost, making it suitable for downstream tasks such as function approximation. In particular, we propose to perform the latter via optimally-weighted least-squares polynomial approximation in the different active subspaces, and we present an adaptive algorithm to choose dynamically the dimensions of the active subspaces and polynomial spaces. We demonstrate the practical viability of the MLAS method with polynomial approximation through numerical experiments based on random partial differential equations (PDEs).

Nested integration problems arise in various scientific and engineering applications, including Bayesian experimental design, financial risk assessment, and uncertainty quantification. These nested integrals take the form $\int f\left(\int g(\bs{y},\bs{x})\di{}\bs{x}\right)\di{}\bs{y}$, for nonlinear $f$, making them computationally challenging, particularly in high-dimensional settings. Although widely used for single integrals, traditional Monte Carlo (MC) methods can be inefficient when encountering complexities of nested integration. This work introduces a novel multilevel estimator, combining deterministic and randomized quasi-MC (rQMC) methods to handle nested integration problems efficiently. In this context, the inner number of samples and the discretization accuracy of the inner integrand evaluation constitute the level. We provide a comprehensive theoretical analysis of the estimator, deriving error bounds demonstrating significant reductions in bias and variance compared with standard methods. The proposed estimator is particularly effective in scenarios where the integrand is evaluated approximately, as it adapts to different levels of resolution without compromising precision. We verify the performance of our method via numerical experiments, focusing on estimating the expected information gain of experiments. When applied to Gaussian noise in the experiment, a truncation scheme ensures finite error bounds at the same computational complexity as in the bounded noise case up to multiplicative logarithmic terms. The results reveal that the proposed multilevel rQMC estimator outperforms existing MC and rQMC approaches, offering a substantial reduction in computational costs and offering a powerful tool for practitioners dealing with complex, nested integration problems across various domains., Comment: 36 pages, 5 figures

In this article we consider likelihood-based estimation of static parameters for a class of partially observed McKean-Vlasov (POMV) diffusion process with discrete-time observations over a fixed time interval. In particular, using the framework of [5] we develop a new randomized multilevel Monte Carlo method for estimating the parameters, based upon Markovian stochastic approximation methodology. New Markov chain Monte Carlo algorithms for the POMV model are introduced facilitating the application of [5]. We prove, under assumptions, that the expectation of our estimator is biased, but with expected small and controllable bias. Our approach is implemented on several examples.

This paper presents an enhanced adaptive random Fourier features (ARFF) training algorithm for shallow neural networks, building upon the work introduced in "Adaptive Random Fourier Features with Metropolis Sampling", Kammonen et al., \emph{Foundations of Data Science}, 2(3):309--332, 2020. This improved method uses a particle filter-type resampling technique to stabilize the training process and reduce the sensitivity to parameter choices. The Metropolis test can also be omitted when resampling is used, reducing the number of hyperparameters by one and reducing the computational cost per iteration compared to the ARFF method. We present comprehensive numerical experiments demonstrating the efficacy of the proposed algorithm in function regression tasks as a stand-alone method and as a pretraining step before gradient-based optimization, using the Adam optimizer. Furthermore, we apply the proposed algorithm to a simple image regression problem, illustrating its utility in sampling frequencies for the random Fourier features (RFF) layer of coordinate-based multilayer perceptrons. In this context, we use the proposed algorithm to sample the parameters of the RFF layer in an automated manner., Comment: 41 pages

This work develops a particle system addressing the approximation of McKean-Vlasov stochastic differential equations (SDEs). The novelty of the approach lies in involving low discrepancy sequences nontrivially in the construction of a particle system with coupled noise and initial conditions. Weak convergence for SDEs with additive noise is proven. A numerical study demonstrates that the novel approach presented here doubles the respective convergence rates for weak and strong approximation of the mean-field limit, compared with the standard particle system. These rates are proven in the simplified setting of a mean-field ordinary differential equation in terms of appropriate bounds involving the star discrepancy for low discrepancy sequences with a group structure, such as Rank-1 lattice points. This construction nontrivially provides an antithetic multilevel quasi-Monte Carlo estimator. An asymptotic error analysis reveals that the proposed approach outperforms methods based on the classic particle system with independent initial conditions and noise.

This work addresses stochastic optimal control problems where the unknown state evolves in continuous time while partial, noisy, and possibly controllable measurements are only available in discrete time. We develop a framework for controlling such systems, focusing on the measure-valued process of the system's state and the control actions that depend on noisy and incomplete data. Our approach uses a stochastic optimal control framework with a probability measure-valued state, which accommodates noisy measurements and integrates them into control decisions through a Bayesian update mechanism. We characterize the control optimality in terms of a sequence of interlaced Hamilton Jacobi Bellman (HJB) equations coupled with controlled impulse steps at the measurement times. For the case of Gaussian-controlled processes, we derive an equivalent HJB equation whose state variable is finite-dimensional, namely the state's mean and covariance. We demonstrate the effectiveness of our methods through numerical examples. These include control under perfect observations, control under no observations, and control under noisy observations. Our numerical results highlight significant differences in the control strategies and their performance, emphasizing the challenges and computational demands of dealing with uncertainty in state observation.

In this work we employ importance sampling (IS) techniques to track a small over-threshold probability of a running maximum associated with the solution of a stochastic differential equation (SDE) within the framework of ensemble Kalman filtering (EnKF). Between two observation times of the EnKF, we propose to use IS with respect to the initial condition of the SDE, IS with respect to the Wiener process via a stochastic optimal control formulation, and combined IS with respect to both initial condition and Wiener process. Both IS strategies require the approximation of the solution of Kolmogorov Backward equation (KBE) with boundary conditions. In multidimensional settings, we employ a Markovian projection dimension reduction technique to obtain an approximation of the solution of the KBE by just solving a one dimensional PDE. The proposed ideas are tested on three illustrative examples: Double Well SDE, Langevin dynamics and noisy Charney-deVore model, and showcase a significant variance reduction compared to the standard Monte Carlo method and another sampling-based IS technique, namely, multilevel cross entropy., Comment: 35 pages, 13 figures, 5 tables

Efficiently pricing multi-asset options poses a significant challenge in quantitative finance. Fourier methods leverage the regularity properties of the integrand in the Fourier domain to accurately and rapidly value options that typically lack regularity in the physical domain. However, most of the existing Fourier approaches face hurdles in high-dimensional settings due to the tensor product (TP) structure of the commonly employed numerical quadrature techniques. To overcome this difficulty, this work advocates using the randomized quasi-MC (RQMC) quadrature to improve the scalability of Fourier methods with high dimensions. The RQMC technique benefits from the smoothness of the integrand and alleviates the curse of dimensionality while providing practical error estimates. Nonetheless, the applicability of RQMC on the unbounded domain, $\mathbb{R}^d$, requires a domain transformation to $[0,1]^d$, which may result in singularities of the transformed integrand at the corners of the hypercube, and hence deteriorate the performance of RQMC. To circumvent this difficulty, we design an efficient domain transformation procedure based on boundary growth conditions on the transformed integrand. The proposed transformation preserves sufficient regularity of the original integrand for fast convergence of the RQMC method. To validate our analysis, we demonstrate the efficiency of employing RQMC with an appropriate transformation to evaluate options in the Fourier space for various pricing models, payoffs, and dimensions. Finally, we highlight the computational advantage of applying RQMC over MC or TP in the Fourier domain, and over MC in the physical domain for options with up to 15 assets.

We investigate the applicability of the well-known multilevel Monte Carlo (MLMC) method to the class of density-driven flow problems, in particular the problem of salinisation of coastal aquifers. As a test case, we solve the uncertain Henry saltwater intrusion problem. Unknown porosity, permeability and recharge parameters are modelled by using random fields. The classical deterministic Henry problem is non-linear and time-dependent, and can easily take several hours of computing time. Uncertain settings require the solution of multiple realisations of the deterministic problem, and the total computational cost increases drastically. Instead of computing of hundreds random realisations, typically the mean value and the variance are computed. The standard methods such as the Monte Carlo or surrogate-based methods is a good choice, but they compute all stochastic realisations on the same, often, very fine mesh. They also do not balance the stochastic and discretisation errors. These facts motivated us to apply the MLMC method. We demonstrate that by solving the Henry problem on multi-level spatial and temporal meshes, the MLMC method reduces the overall computational and storage costs. To reduce the computing cost further, parallelization is performed in both physical and stochastic spaces. To solve each deterministic scenario, we run the parallel multigrid solver ug4 in a black-box fashion., Comment: 20 pages, 11 figures, 3 tables. arXiv admin note: substantial text overlap with arXiv:2302.07804

This paper addresses the difficulty of characterizing the time-varying nature of fading channels. The current time-invariant models often fall short of capturing and tracking these dynamic characteristics. To overcome this limitation, we explore using of stochastic differential equations (SDEs) and Markovian projection to model signal envelope variations, considering scenarios involving Rayleigh, Rice, and Hoyt distributions. Furthermore, it is of practical interest to study the performance of channels modeled by SDEs. In this work, we investigate the fade duration metric, representing the time during which the signal remains below a specified threshold within a fixed time interval. We estimate the complementary cumulative distribution function (CCDF) of the fade duration using Monte Carlo simulations, and analyze the influence of system parameters on its behavior. Finally, we leverage importance sampling, a known variance-reduction technique, to estimate the tail of the CCDF efficiently.

We present experimental results highlighting two key differences resulting from the choice of training algorithm for two-layer neural networks. The spectral bias of neural networks is well known, while the spectral bias dependence on the choice of training algorithm is less studied. Our experiments demonstrate that an adaptive random Fourier features algorithm (ARFF) can yield a spectral bias closer to zero compared to the stochastic gradient descent optimizer (SGD). Additionally, we train two identically structured classifiers, employing SGD and ARFF, to the same accuracy levels and empirically assess their robustness against adversarial noise attacks., Comment: 6 Pages, 4 Figures; Accepted in the International Conference on Scientific Computing and Machine Learning

This work considers a short-term, continuous time setting characterized by a coupled power supply system controlled exclusively by a single provider and comprising a cascade of hydropower systems (dams), fossil fuel power stations, and a storage capacity modeled by a single large battery. Cascaded hydropower generators introduce time-delay effects in the state dynamics, which are modeled with differential equations, making it impossible to use classical dynamic programming. We address this issue by introducing a novel Lagrangian relaxation technique over continuous-time constraints, constructing a nearly optimal policy efficiently. This approach yields a convex, nonsmooth optimization dual problem to recover the optimal Lagrangian multipliers, which is numerically solved using a limited memory bundle method. At each step of the dual optimization, we need to solve an optimization subproblem. Given the current values of the Lagrangian multipliers, the time delays are no longer active, and we can solve a corresponding nonlinear Hamilton--Jacobi--Bellman (HJB) Partial Differential Equation (PDE) for the optimization subproblem. The HJB PDE solver provides both the current value of the dual function and its subgradient, and is trivially parallelizable over the state space for each time step. To handle the infinite-dimensional nature of the Lagrange multipliers, we design an adaptive refinement strategy to control the duality gap. Furthermore, we use a penalization technique for the constructed admissible primal solution to smooth the controls while achieving a sufficiently small duality gap. Numerical results based on the Uruguayan power system demonstrate the efficiency of the proposed mathematical models and numerical approach.

This study analyzes the nonasymptotic convergence behavior of the quasi-Monte Carlo (QMC) method with applications to linear elliptic partial differential equations (PDEs) with lognormal coefficients. Building upon the error analysis presented in (Owen, 2006), we derive a nonasymptotic convergence estimate depending on the specific integrands, the input dimensionality, and the finite number of samples used in the QMC quadrature. We discuss the effects of the variance and dimensionality of the input random variable. Then, we apply the QMC method with importance sampling (IS) to approximate deterministic, real-valued, bounded linear functionals that depend on the solution of a linear elliptic PDE with a lognormal diffusivity coefficient in bounded domains of $\mathbb{R}^d$, where the random coefficient is modeled as a stationary Gaussian random field parameterized by the trigonometric and wavelet-type basis. We propose two types of IS distributions, analyze their effects on the QMC convergence rate, and observe the improvements., Comment: 47 pages, 9 figures

Finding the optimal design of experiments in the Bayesian setting typically requires estimation and optimization of the expected information gain functional. This functional consists of one outer and one inner integral, separated by the logarithm function applied to the inner integral. When the mathematical model of the experiment contains uncertainty about the parameters of interest and nuisance uncertainty, (i.e., uncertainty about parameters that affect the model but are not themselves of interest to the experimenter), two inner integrals must be estimated. Thus, the already considerable computational effort required to determine good approximations of the expected information gain is increased further. The Laplace approximation has been applied successfully in the context of experimental design in various ways, and we propose two novel estimators featuring the Laplace approximation to alleviate the computational burden of both inner integrals considerably. The first estimator applies Laplace's method followed by a Laplace approximation, introducing a bias. The second estimator uses two Laplace approximations as importance sampling measures for Monte Carlo approximations of the inner integrals. Both estimators use Monte Carlo approximation for the remaining outer integral estimation. We provide four numerical examples demonstrating the applicability and effectiveness of our proposed estimators., Comment: 20 pages, 10 figures

In this work, we consider the general problem of constructing a neural network surrogate model using multi-fidelity information. Motivated by error-complexity estimates for ReLU neural networks, we formulate the correlation between an inexpensive low-fidelity model and an expensive high-fidelity model as a possibly non-linear residual function. This function defines a mapping between 1) the shared input space of the models along with the low-fidelity model output, and 2) the discrepancy between the outputs of the two models. The computational framework proceeds by training two neural networks to work in concert. The first network learns the residual function on a small set of high- and low-fidelity data. Once trained, this network is used to generate additional synthetic high-fidelity data, which is used in the training of the second network. The trained second network then acts as our surrogate for the high-fidelity quantity of interest. We present four numerical examples to demonstrate the power of the proposed framework, showing that significant savings in computational cost may be achieved when the output predictions are desired to be accurate within small tolerances.

This work introduces a novel approach that combines the multi-index Monte Carlo (MC) method with importance sampling (IS) to estimate rare event quantities expressed as an expectation of a smooth observable of solutions to a broad class of McKean-Vlasov stochastic differential equations. We extend the double loop Monte Carlo (DLMC) estimator, previously introduced in our works (Ben Rached et al., 2022a,b), to the multi-index setting. We formulate a new multi-index DLMC estimator and conduct a comprehensive cost-error analysis, leading to improved complexity results. To address rare events, an importance sampling scheme is applied using stochastic optimal control of the single level DLMC estimator. This combination of IS and multi-index DLMC not only reduces computational complexity by two orders but also significantly decreases the associated constant compared to vanilla MC. The effectiveness of the proposed multi-index DLMC estimator is demonstrated using the Kuramoto model from statistical physics. The results confirm a reduced complexity from $\mathcal{O}(\mathrm{TOL}_{\mathrm{r}}^{-4})$ for the single level DLMC estimator (Ben Rached et al., 2022a) to $\mathcal{O}(\mathrm{TOL}_{\mathrm{r}}^{-2} (\log \mathrm{TOL}_{\mathrm{r}}^{-1})^2)$ for the considered example, while ensuring accurate estimation of rare event quantities within the prescribed relative error tolerance $\mathrm{TOL}_\mathrm{r}$., Comment: Extension to works 2207.06926 and 2208.03225

We address the computational efficiency in solving the A-optimal Bayesian design of experiments problems for which the observational map is based on partial differential equations and, consequently, is computationally expensive to evaluate. A-optimality is a widely used and easy-to-interpret criterion for Bayesian experimental design. This criterion seeks the optimal experimental design by minimizing the expected conditional variance, which is also known as the expected posterior variance. This study presents a novel likelihood-free approach to the A-optimal experimental design that does not require sampling or integrating the Bayesian posterior distribution. The expected conditional variance is obtained via the variance of the conditional expectation using the law of total variance, and we take advantage of the orthogonal projection property to approximate the conditional expectation. We derive an asymptotic error estimation for the proposed estimator of the expected conditional variance and show that the intractability of the posterior distribution does not affect the performance of our approach. We use an artificial neural network (ANN) to approximate the nonlinear conditional expectation in the implementation of our method. We then extend our approach for dealing with the case that the domain of experimental design parameters is continuous by integrating the training process of the ANN into minimizing the expected conditional variance. Through numerical experiments, we demonstrate that our method greatly reduces the number of observation model evaluations compared with widely used importance sampling-based approaches. This reduction is crucial, considering the high computational cost of the observational models. Code is available at https://github.com/vinh-tr-hoang/DOEviaPACE.

We propose a novel alternative approach to our previous work (Ben Hammouda et al., 2023) to improve the efficiency of Monte Carlo (MC) estimators for rare event probabilities for stochastic reaction networks (SRNs). In the same spirit of (Ben Hammouda et al., 2023), an efficient path-dependent measure change is derived based on a connection between determining optimal importance sampling (IS) parameters within a class of probability measures and a stochastic optimal control formulation, corresponding to solving a variance minimization problem. In this work, we propose a novel approach to address the encountered curse of dimensionality by mapping the problem to a significantly lower-dimensional space via a Markovian projection (MP) idea. The output of this model reduction technique is a low-dimensional SRN (potentially even one dimensional) that preserves the marginal distribution of the original high-dimensional SRN system. The dynamics of the projected process are obtained by solving a related optimization problem via a discrete $L^2$ regression. By solving the resulting projected Hamilton-Jacobi-Bellman (HJB) equations for the reduced-dimensional SRN, we obtain projected IS parameters, which are then mapped back to the original full-dimensional SRN system, resulting in an efficient IS-MC estimator for rare events probabilities of the full-dimensional SRN. Our analysis and numerical experiments reveal that the proposed MP-HJB-IS approach substantially reduces the MC estimator variance, resulting in a lower computational complexity in the rare event regime than standard MC estimators.

This work employs the Birnbaum--Saunders distribution to model the fatigue life of metallic materials under cyclic loading and compares it with the normal distribution. Fatigue-limit models are fitted to three datasets of unnotched specimens of 75S-T6 aluminum alloys and carbon laminate with different loading types. A new equivalent stress definition that accounts for the effect of the experiment type is proposed. The results show that the Birnbaum--Saunders distribution consistently outperforms the normal distribution in fitting the fatigue data and provides more accurate predictions of fatigue life and survival probability.

Nested integration of the form $\int f\left(\int g(\bs{y},\bs{x})\di{}\bs{x}\right)\di{}\bs{y}$, characterized by an outer integral connected to an inner integral through a nonlinear function $f$, is a challenging problem in various fields, such as engineering and mathematical finance. The available numerical methods for nested integration based on Monte Carlo (MC) methods can be prohibitively expensive owing to the error propagating from the inner to the outer integral. Attempts to enhance the efficiency of these approximations using the quasi-MC (QMC) or randomized QMC (rQMC) method have focused on either the inner or outer integral approximation. This work introduces a novel nested rQMC method that simultaneously addresses the approximation of the inner and outer integrals. The method leverages the unique nested integral structure to offer a more efficient approximation mechanism. Incorporating Owen's scrambling techniques, we address integrands exhibiting infinite variation in the Hardy--Krause sense, enabling theoretically sound error estimates. As the primary contribution, we derive asymptotic error bounds for the bias and variance of our estimator, along with the regularity conditions under which these bounds can be attained. Moreover, we derive a truncation scheme for applications in the context of expected information gain estimation and indicate how to use importance sampling to remedy the measure concentration arising in the inner integral. We verify the estimator quality through numerical experiments by comparing the computational efficiency of the nested rQMC method against standard nested MC estimation for two case studies: one in thermomechanics and the other in pharmacokinetics. These examples highlight the computational savings and enhanced applicability of the proposed approach., Comment: 29 pages, 9 figures

In this work, we extend the data-driven It\^{o} stochastic differential equation (SDE) framework for the pathwise assessment of short-term forecast errors to account for the time-dependent upper bound that naturally constrains the observable historical data and forecast. We propose a new nonlinear and time-inhomogeneous SDE model with a Jacobi-type diffusion term for the phenomenon of interest, simultaneously driven by the forecast and the constraining upper bound. We rigorously demonstrate the existence and uniqueness of a strong solution to the SDE model by imposing a condition for the time-varying mean-reversion parameter appearing in the drift term. The normalized forecast function is thresholded to keep such mean-reversion parameters bounded. The SDE model parameter calibration is applied to user-selected approximations of the likelihood function. Another novel contribution is estimating the unknown transition density of the forecast error process with a tailored kernel smoothing technique with the control variate method, coupling an adequate SDE to the original one. We provide a theoretical study about how to choose the optimal bandwidth. We fit the model to the 2019 photovoltaic (PV) solar power daily production and forecast data in Uruguay, computing the daily maximum solar PV production estimation. Two statistical versions of the constrained SDE model are fit, with the beta and truncated normal distributions as proxies for the transition density. Empirical results include simulations of the normalized solar PV power production and pathwise confidence bands generated through an indirect inference method. An objective comparison of optimal parametric points associated with the two selected statistical approximations is provided by applying our innovative kernel smoothing estimation technique of the transition function of the forecast error process., Comment: 28 pages, 18 figures

In this work, we further develop the Physics-informed Spectral Learning (PiSL) by Espath et al. \cite{Esp21} based on a discrete $L^2$ projection to solve the discrete Hodge--Helmholtz decomposition from sparse data. Within this physics-informed statistical learning framework, we adaptively build a sparse set of Fourier basis functions with corresponding coefficients by solving a sequence of minimization problems where the set of basis functions is augmented greedily at each optimization problem. Moreover, our PiSL computational framework enjoys spectral (exponential) convergence. We regularize the minimization problems with the seminorm of the fractional Sobolev space in a Tikhonov fashion. In the Fourier setting, the divergence- and curl-free constraints become a finite set of linear algebraic equations. The proposed computational framework combines supervised and unsupervised learning techniques in that we use data concomitantly with the projection onto divergence- and curl-free spaces. We assess the capabilities of our method in various numerical examples including the `Storm of the Century' with satellite data from 1993.

We consider a class of density-driven flow problems. We are particularly interested in the problem of the salinization of coastal aquifers. We consider the Henry saltwater intrusion problem with uncertain porosity, permeability, and recharge parameters as a test case. The reason for the presence of uncertainties is the lack of knowledge, inaccurate measurements, and inability to measure parameters at each spatial or time location. This problem is nonlinear and time-dependent. The solution is the salt mass fraction, which is uncertain and changes in time. Uncertainties in porosity, permeability, recharge, and mass fraction are modeled using random fields. This work investigates the applicability of the well-known multilevel Monte Carlo (MLMC) method for such problems. The MLMC method can reduce the total computational and storage costs. Moreover, the MLMC method runs multiple scenarios on different spatial and time meshes and then estimates the mean value of the mass fraction. The parallelization is performed in both the physical space and stochastic space. To solve every deterministic scenario, we run the parallel multigrid solver ug4 in a black-box fashion. We use the solution obtained from the quasi-Monte Carlo method as a reference solution., Comment: 24 pages, 3 tables, 11 figures

In this study, we propose a new numerical scheme for physics-informed neural networks (PINNs) that enables precise and inexpensive solution for partial differential equations (PDEs) in case of arbitrary geometries while strictly enforcing Dirichlet boundary conditions. The proposed approach combines admissible NURBS parametrizations required to define the physical domain and the Dirichlet boundary conditions with a PINN solver. The fundamental boundary conditions are automatically satisfied in this novel Deep NURBS framework. We verified our new approach using two-dimensional elliptic PDEs when considering arbitrary geometries, including non-Lipschitz domains. Compared to the classical PINN solver, the Deep NURBS estimator has a remarkably high convergence rate for all the studied problems. Moreover, a desirable accuracy was realized for most of the studied PDEs using only one hidden layer of neural networks. This novel approach is considered to pave the way for more effective solutions for high-dimensional problems by allowing for more realistic physics-informed statistical learning to solve PDE-based variational problems.

This work combines multilevel Monte Carlo (MLMC) with importance sampling to estimate rare-event quantities that can be expressed as the expectation of a Lipschitz observable of the solution to a broad class of McKean--Vlasov stochastic differential equations. We extend the double loop Monte Carlo (DLMC) estimator introduced in this context in (Ben Rached et al., 2023) to the multilevel setting. We formulate a novel multilevel DLMC estimator and perform a comprehensive cost-error analysis yielding new and improved complexity results. Crucially, we devise an antithetic sampler to estimate level differences guaranteeing reduced computational complexity for the multilevel DLMC estimator compared with the single-level DLMC estimator. To address rare events, we apply the importance sampling scheme, obtained via stochastic optimal control in (Ben Rached et al., 2023), over all levels of the multilevel DLMC estimator. Combining importance sampling and multilevel DLMC reduces computational complexity by one order and drastically reduces the associated constant compared to the single-level DLMC estimator without importance sampling. We illustrate the effectiveness of the proposed multilevel DLMC estimator on the Kuramoto model from statistical physics with Lipschitz observables, confirming the reduced complexity from $\mathcal{O}(\mathrm{TOL}_{\mathrm{r}}^{-4})$ for the single-level DLMC estimator to $\mathcal{O}(\mathrm{TOL}_{\mathrm{r}}^{-3})$ while providing a feasible estimate of rare-event quantities up to prescribed relative error tolerance $\mathrm{TOL}_{\mathrm{r}}$., Comment: Follow-up to arXiv:2207.06926

Using quasi-Newton methods in stochastic optimization is not a trivial task given the difficulty of extracting curvature information from the noisy gradients. Moreover, pre-conditioning noisy gradient observations tend to amplify the noise. We propose a Bayesian approach to obtain a Hessian matrix approximation for stochastic optimization that minimizes the secant equations residue while retaining the extreme eigenvalues between a specified range. Thus, the proposed approach assists stochastic gradient descent to converge to local minima without augmenting gradient noise. We propose maximizing the log posterior using the Newton-CG method. Numerical results on a stochastic quadratic function and an $\ell_2$-regularized logistic regression problem are presented. In all the cases tested, our approach improves the convergence of stochastic gradient descent, compensating for the overhead of solving the log posterior maximization. In particular, pre-conditioning the stochastic gradient with the inverse of our Hessian approximation becomes more advantageous the larger the condition number of the problem is.

This paper investigates Monte Carlo (MC) methods to estimate probabilities of rare events associated with solutions to the $d$-dimensional McKean-Vlasov stochastic differential equation (MV-SDE). MV-SDEs are usually approximated using a stochastic interacting $P$-particle system, which is a set of $P$ coupled $d$-dimensional stochastic differential equations (SDEs). Importance sampling (IS) is a common technique for reducing high relative variance of MC estimators of rare-event probabilities. We first derive a zero-variance IS change of measure for the quantity of interest by using stochastic optimal control theory. However, when this change of measure is applied to stochastic particle systems, it yields a $P \times d$-dimensional partial differential control equation (PDE), which is computationally expensive to solve. To address this issue, we use the decoupling approach introduced in [dos Reis et al., 2023], generating a $d$-dimensional control PDE for a zero-variance estimator of the decoupled SDE. Based on this approach, we develop a computationally efficient double loop MC (DLMC) estimator. We conduct a comprehensive numerical error and work analysis of the DLMC estimator. As a result, we show optimal complexity of $\mathcal{O}(\mathrm{TOL}_{\mathrm{r}}^{-4})$ with a significantly reduced constant to achieve a prescribed relative error tolerance $\mathrm{TOL}_{\mathrm{r}}$. Subsequently, we propose an adaptive DLMC method combined with IS to numerically estimate rare-event probabilities, substantially reducing relative variance and computational runtimes required to achieve a given $\mathrm{TOL}_{\mathrm{r}}$ compared with standard MC estimators in the absence of IS. Numerical experiments are performed on the Kuramoto model from statistical physics.

We present an adaptive multilevel Monte Carlo (AMLMC) algorithm for approximating deterministic, real-valued, bounded linear functionals that depend on the solution of a linear elliptic PDE with a lognormal diffusivity coefficient and geometric singularities in bounded domains of $\mathbb{R}^d$. Our AMLMC algorithm is built on the results of the weak convergence rates in the work [Moon et al., BIT Numer. Math., 46 (2006), 367-407] for an adaptive algorithm using isoparametric d-linear quadrilateral finite element approximations and the dual weighted residual error representation in a deterministic setting. Designed to suit the geometric nature of the singularities in the solution, our AMLMC algorithm uses a sequence of deterministic, non-uniform auxiliary meshes as a building block. The deterministic adaptive algorithm generates these meshes, corresponding to a geometrically decreasing sequence of tolerances. For a given realization of the diffusivity coefficient and accuracy level, AMLMC constructs its approximate sample using the first mesh in the hierarchy that satisfies the corresponding bias accuracy constraint. This adaptive approach is particularly useful for the lognormal case treated here, which lacks uniform coercivity and thus produces functional outputs that vary over orders of magnitude when sampled. We discuss iterative solvers and compare their efficiency with direct ones. To reduce computational work, we propose a stopping criterion for the iterative solver with respect to the quantity of interest, the realization of the diffusivity coefficient, and the desired level of AMLMC approximation. From the numerical experiments, based on a Fourier expansion of the coefficient field, we observe improvements in efficiency compared with both standard Monte Carlo and standard MLMC for a problem with a singularity similar to that at the tip of a slit modeling a crack., Comment: Fixed typos, included more proofs and discussions

Efficiently pricing multi-asset options is a challenging problem in quantitative finance. When the characteristic function is available, Fourier-based methods are competitive compared to alternative techniques because the integrand in the frequency space often has a higher regularity than that in the physical space. However, when designing a numerical quadrature method for most Fourier pricing approaches, two key aspects affecting the numerical complexity should be carefully considered: (i) the choice of damping parameters that ensure integrability and control the regularity class of the integrand and (ii) the effective treatment of high dimensionality. We propose an efficient numerical method for pricing European multi-asset options based on two complementary ideas to address these challenges. First, we smooth the Fourier integrand via an optimized choice of the damping parameters based on a proposed optimization rule. Second, we employ sparsification and dimension-adaptivity techniques to accelerate the convergence of the quadrature in high dimensions. The extensive numerical study on basket and rainbow options under the multivariate geometric Brownian motion and some L\'evy models demonstrates the advantages of adaptivity and the damping rule on the numerical complexity of quadrature methods. Moreover, for the tested two-asset examples, the proposed approach outperforms the COS method in terms of computational time. Finally, we show significant speed-up compared to the Monte Carlo method for up to six dimensions.

To model manifold data using normalizing flows, we employ isometric autoencoders to design embeddings with explicit inverses that do not distort the probability distribution. Using isometries separates manifold learning and density estimation and enables training of both parts to high accuracy. Thus, model selection and tuning are simplified compared to existing injective normalizing flows. Applied to data sets on (approximately) flat manifolds, the combined approach generates high-quality data., Comment: 11 pages, 7 figures, 4 tables

Estimating the expectations of functionals applied to sums of random variables (RVs) is a well-known problem encountered in many challenging applications. Generally, closed-form expressions of these quantities are out of reach. A naive Monte Carlo simulation is an alternative approach. However, this method requires numerous samples for rare event problems. Therefore, it is paramount to use variance reduction techniques to develop fast and efficient estimation methods. In this work, we use importance sampling (IS), known for its efficiency in requiring fewer computations to achieve the same accuracy requirements. We propose a state-dependent IS scheme based on a stochastic optimal control formulation, where the control is dependent on state and time. We aim to calculate rare event quantities that could be written as an expectation of a functional of the sums of independent RVs. The proposed algorithm is generic and can be applied without restrictions on the univariate distributions of RVs or the functional applied to the sum. We apply this approach to the log-normal distribution to compute the left tail and cumulative distribution of the ratio of independent RVs. For each case, we numerically demonstrate that the proposed state-dependent IS algorithm compares favorably to most well-known estimators dealing with similar problems.

Calculating the expected information gain in optimal Bayesian experimental design typically relies on nested Monte Carlo sampling. When the model also contains nuisance parameters, which are parameters that contribute to the overall uncertainty of the system but are of no interest in the Bayesian design framework, this introduces a second inner loop. We propose and derive a small-noise approximation for this additional inner loop. The computational cost of our method can be further reduced by applying a Laplace approximation to the remaining inner loop. Thus, we present two methods, the small-noise Double-loop Monte Carlo and small-noise Monte Carlo Laplace methods. Moreover, we demonstrate that the total complexity of these two approaches remains comparable to the case without nuisance uncertainty. To assess the efficiency of these methods, we present three examples, and the last example includes the partial differential equation for the electrical impedance tomography experiment for composite laminate materials., Comment: 22 pages, 11 figures, Added references, Addressed referee suggestions, Corrected equations in Section 5 and Appendix B; No significant impact on numerical results, Updated figures, Corrected notation in Section 6, Corrected typos

When approximating the expectations of a functional of a solution to a stochastic differential equation, the numerical performance of deterministic quadrature methods, such as sparse grid quadrature and quasi-Monte Carlo (QMC) methods, may critically depend on the regularity of the integrand. To overcome this issue and improve the regularity structure of the problem, we consider cases in which analytic smoothing (bias-free mollification) cannot be performed and introduce a novel numerical smoothing approach by combining a root-finding method with a one-dimensional numerical integration with respect to a single well-chosen variable. We prove that, under appropriate conditions, the resulting function of the remaining variables is highly smooth, potentially affording the improved efficiency of adaptive sparse grid quadrature (ASGQ) and QMC methods, particularly when combined with hierarchical transformations (ie., the Brownian bridge and Richardson extrapolation on the weak error). This approach facilitates the effective treatment of high dimensionality. Our study is motivated by option pricing problems, focusing on dynamics where the discretization of the asset price is necessary. Based on our analysis and numerical experiments, we demonstrate the advantages of combining numerical smoothing with the ASGQ and QMC methods over these methods without smoothing and the Monte Carlo approach. Finally, our approach is generic and can be applied to solve a broad class of problems, particularly approximating distribution functions, computing financial Greeks, and estimating risk quantities., Comment: arXiv admin note: substantial text overlap with arXiv:2003.05708

We explore efficient estimation of statistical quantities, particularly rare event probabilities, for stochastic reaction networks. Consequently, we propose an importance sampling (IS) approach to improve the Monte Carlo (MC) estimator efficiency based on an approximate tau-leap scheme. The crucial step in the IS framework is choosing an appropriate change of probability measure to achieve substantial variance reduction. This task is typically challenging and often requires insights into the underlying problem. Therefore, we propose an automated approach to obtain a highly efficient path-dependent measure change based on an original connection in the stochastic reaction network context between finding optimal IS parameters within a class of probability measures and a stochastic optimal control formulation. Optimal IS parameters are obtained by solving a variance minimization problem. First, we derive an associated dynamic programming equation. Analytically solving this backward equation is challenging, hence we propose an approximate dynamic programming formulation to find near-optimal control parameters. To mitigate the curse of dimensionality, we propose a learning-based method to approximate the value function using a neural network, where the parameters are determined via a stochastic optimization algorithm. Our analysis and numerical experiments verify that the proposed learning-based IS approach substantially reduces MC estimator variance, resulting in a lower computational complexity in the rare event regime, compared with standard tau-leap MC estimators.

In this study, we demonstrate that the norm test and inner product/orthogonality test presented in \cite{Bol18} are equivalent in terms of the convergence rates associated with Stochastic Gradient Descent (SGD) methods if $\epsilon^2=\theta^2+\nu^2$ with specific choices of $\theta$ and $\nu$. Here, $\epsilon$ controls the relative statistical error of the norm of the gradient while $\theta$ and $\nu$ control the relative statistical error of the gradient in the direction of the gradient and in the direction orthogonal to the gradient, respectively. Furthermore, we demonstrate that the inner product/orthogonality test can be as inexpensive as the norm test in the best case scenario if $\theta$ and $\nu$ are optimally selected, but the inner product/orthogonality test will never be more computationally affordable than the norm test if $\epsilon^2=\theta^2+\nu^2$. Finally, we present two stochastic optimization problems to illustrate our results.

We reconstruct the velocity field of incompressible flows given a finite set of measurements. For the spatial approximation, we introduce the Sparse Fourier divergence-free (SFdf) approximation based on a discrete $L^2$ projection. Within this physics-informed type of statistical learning framework, we adaptively build a sparse set of Fourier basis functions with corresponding coefficients by solving a sequence of minimization problems where the set of basis functions is augmented greedily at each optimization problem. We regularize our minimization problems with the seminorm of the fractional Sobolev space in a Tikhonov fashion. In the Fourier setting, the incompressibility (divergence-free) constraint becomes a finite set of linear algebraic equations. We couple our spatial approximation with the truncated Singular Value Decomposition (SVD) of the flow measurements for temporal compression. Our computational framework thus combines supervised and unsupervised learning techniques. We assess the capabilities of our method in various numerical examples arising in fluid mechanics.

This paper presents the machine learning-based ensemble conditional mean filter (ML-EnCMF) -- a filtering method based on the conditional mean filter (CMF) previously introduced in the literature. The updated mean of the CMF matches that of the posterior, obtained by applying Bayes' rule on the filter's forecast distribution. Moreover, we show that the CMF's updated covariance coincides with the expected conditional covariance. Implementing the EnCMF requires computing the conditional mean (CM). A likelihood-based estimator is prone to significant errors for small ensemble sizes, causing the filter divergence. We develop a systematical methodology for integrating machine learning into the EnCMF based on the CM's orthogonal projection property. First, we use a combination of an artificial neural network (ANN) and a linear function, obtained based on the ensemble Kalman filter (EnKF), to approximate the CM, enabling the ML-EnCMF to inherit EnKF's advantages. Secondly, we apply a suitable variance reduction technique to reduce statistical errors when estimating loss function. Lastly, we propose a model selection procedure for element-wisely selecting the applied filter, i.e., either the EnKF or ML-EnCMF, at each updating step. We demonstrate the ML-EnCMF performance using the Lorenz-63 and Lorenz-96 systems and show that the ML-EnCMF outperforms the EnKF and the likelihood-based EnCMF.

In this work, we present, analyze, and implement a class of Multi-Level Markov chain Monte Carlo (ML-MCMC) algorithms based on independent Metropolis-Hastings proposals for Bayesian inverse problems. In this context, the likelihood function involves solving a complex differential model, which is then approximated on a sequence of increasingly accurate discretizations. The key point of this algorithm is to construct highly coupled Markov chains together with the standard Multi-level Monte Carlo argument to obtain a better cost-tolerance complexity than a single-level MCMC algorithm. Our method extends the ideas of Dodwell, et al. "A hierarchical multilevel Markov chain Monte Carlo algorithm with applications to uncertainty quantification in subsurface flow," \textit{SIAM/ASA Journal on Uncertainty Quantification 3.1 (2015): 1075-1108,} to a wider range of proposal distributions. We present a thorough convergence analysis of the ML-MCMC method proposed, and show, in particular, that (i) under some mild conditions on the (independent) proposals and the family of posteriors, there exists a unique invariant probability measure for the coupled chains generated by our method, and (ii) that such coupled chains are uniformly ergodic. We also generalize the cost-tolerance theorem of Dodwell et al., to our wider class of ML-MCMC algorithms. Finally, we propose a self-tuning continuation-type ML-MCMC algorithm (C-ML-MCMC). The presented method is tested on an array of academic examples, where some of our theoretical results are numerically verified. These numerical experiments evidence how our extended ML-MCMC method is robust when targeting some \emph{pathological} posteriors, for which some of the previously proposed ML-MCMC algorithms fail.

Neural networks-based learning of the distribution of non-dispatchable renewable electricity generation from sources such as photovoltaics (PV) and wind as well as load demands has recently gained attention. Normalizing flow density models are particularly well suited for this task due to the training through direct log-likelihood maximization. However, research from the field of image generation has shown that standard normalizing flows can only learn smeared-out versions of manifold distributions. Previous works on normalizing flow-based scenario generation do not address this issue, and the smeared-out distributions result in the sampling of noisy time series. In this paper, we exploit the isometry of the principal component analysis (PCA), which sets up the normalizing flow in a lower-dimensional space while maintaining the direct and computationally efficient likelihood maximization. We train the resulting principal component flow (PCF) on data of PV and wind power generation as well as load demand in Germany in the years 2013 to 2015. The results of this investigation show that the PCF preserves critical features of the original distributions, such as the probability density and frequency behavior of the time series. The application of the PCF is, however, not limited to renewable power generation but rather extends to any data set, time series, or otherwise, which can be efficiently reduced using PCA., Comment: 18 pages, 7 figures

In this work we combine ideas from multi-index Monte Carlo and ensemble Kalman filtering (EnKF) to produce a highly efficient filtering method called multi-index EnKF (MIEnKF). MIEnKF is based on independent samples of four-coupled EnKF estimators on a multi-index hierarchy of resolution levels, and it may be viewed as an extension of the multilevel EnKF (MLEnKF) method developed by the same authors in 2020. Multi-index here refers to a two-index method, consisting of a hierarchy of EnKF estimators that are coupled in two degrees of freedom: time discretization and ensemble size. Under certain assumptions, when strong coupling between solutions on neighboring numerical resolutions is attainable, the MIEnKF method is proven to be more tractable than EnKF and MLEnKF. Said efficiency gains are also verified numerically in a series of test problems., Comment: 29 pages, 11 figures, 1 algorithm. Accepted by the Journal of Computational Physics. The computer code can be downloaded from https://github.com/GaukharSH/mienkf

We investigate the use of spatial interpolation methods for reconstructing the horizontal near-surface wind field given a sparse set of measurements. In particular, random Fourier features is compared to a set of benchmark methods including Kriging and Inverse distance weighting. Random Fourier features is a linear model $\beta(\pmb x) = \sum_{k=1}^K \beta_k e^{i\omega_k \pmb x}$ approximating the velocity field, with frequencies $\omega_k$ randomly sampled and amplitudes $\beta_k$ trained to minimize a loss function. We include a physically motivated divergence penalty term $|\nabla \cdot \beta(\pmb x)|^2$, as well as a penalty on the Sobolev norm. We derive a bound on the generalization error and derive a sampling density that minimizes the bound. Following (arXiv:2007.10683 [math.NA]), we devise an adaptive Metropolis-Hastings algorithm for sampling the frequencies of the optimal distribution. In our experiments, our random Fourier features model outperforms the benchmark models.

We discuss estimating the probability that the sum of nonnegative independent and identically distributed random variables falls below a given threshold, i.e., $\mathbb{P}(\sum_{i=1}^{N}{X_i} \leq \gamma)$, via importance sampling (IS). We are particularly interested in the rare event regime when $N$ is large and/or $\gamma$ is small. The exponential twisting is a popular technique for similar problems that, in most cases, compares favorably to other estimators. However, it has some limitations: i) it assumes the knowledge of the moment generating function of $X_i$ and ii) sampling under the new IS PDF is not straightforward and might be expensive. The aim of this work is to propose an alternative IS PDF that approximately yields, for certain classes of distributions and in the rare event regime, at least the same performance as the exponential twisting technique and, at the same time, does not introduce serious limitations. The first class includes distributions whose probability density functions (PDFs) are asymptotically equivalent, as $x \rightarrow 0$, to $bx^{p}$, for $p>-1$ and $b>0$. For this class of distributions, the Gamma IS PDF with appropriately chosen parameters retrieves approximately, in the rare event regime corresponding to small values of $\gamma$ and/or large values of $N$, the same performance of the estimator based on the use of the exponential twisting technique. In the second class, we consider the Log-normal setting, whose PDF at zero vanishes faster than any polynomial, and we show numerically that a Gamma IS PDF with optimized parameters clearly outperforms the exponential twisting IS PDF. Numerical experiments validate the efficiency of the proposed estimator in delivering a highly accurate estimate in the regime of large $N$ and/or small $\gamma$.

We present a proof showing that the weak error of a system of $n$ interacting stochastic particles approximating the solution of the McKean-Vlasov equation is $\mathcal O(n^{-1})$. Our proof is based on the Kolmogorov backward equation for the particle system and bounds on the derivatives of its solution, which we derive more generally using the variations of the stochastic particle system. The convergence rate is verified by numerical experiments, which also indicate that the assumptions made here and in the literature can be relaxed.

We here introduce Multi-Iteration Stochastic Optimizers, a novel class of first-order stochastic optimizers where the relative $L^2$ error is estimated and controlled using successive control variates along the path of iterations. By exploiting the correlation between iterates, control variates may reduce the estimator's variance so that an accurate estimation of the mean gradient becomes computationally affordable. We name the estimator of the mean gradient Multi-Iteration stochastiC Estimator (MICE). In principle, MICE can be flexibly coupled with any first-order stochastic optimizer, given its non-intrusive nature. Our generic algorithm adaptively decides which iterates to keep in its index set. We present an error analysis of MICE and a convergence analysis of Multi-Iteration Stochastic Optimizers for different classes of problems, including some non-convex cases. Within the smooth, strongly convex setting, we show that to approximate a minimizer with accuracy $tol$, SGD-MICE requires, on average, $O(tol^{-1})$ stochastic gradient evaluations, while SGD with adaptive batch sizes requires $O(tol^{-1} \log(tol^{-1}))$, correspondingly. Moreover, in a numerical evaluation, SGD-MICE achieved tol with less than 3% the number of gradient evaluations than adaptive batch SGD. The MICE estimator provides a straightforward stopping criterion based on the gradient norm that is validated in consistency tests. To assess the efficiency of MICE, we present several examples in which we use SGD-MICE and Adam-MICE. We include one example based on a stochastic adaptation of the Rosenbrock function and logistic regression training for various datasets. When compared to SGD, SAG, SAGA, SVRG, and SARAH, the Multi-Iteration Stochastic Optimizers reduced, without the need to tune parameters for each example, the gradient sampling cost in all cases tested, also being competitive in runtime in some cases.

Estimates of the generalization error are proved for a residual neural network with $L$ random Fourier features layers $\bar z_{\ell+1}=\bar z_\ell + \mathrm{Re}\sum_{k=1}^K\bar b_{\ell k}e^{\mathrm{i}\omega_{\ell k}\bar z_\ell}+ \mathrm{Re}\sum_{k=1}^K\bar c_{\ell k}e^{\mathrm{i}\omega'_{\ell k}\cdot x}$. An optimal distribution for the frequencies $(\omega_{\ell k},\omega'_{\ell k})$ of the random Fourier features $e^{\mathrm{i}\omega_{\ell k}\bar z_\ell}$ and $e^{\mathrm{i}\omega'_{\ell k}\cdot x}$ is derived. This derivation is based on the corresponding generalization error for the approximation of the function values $f(x)$. The generalization error turns out to be smaller than the estimate ${\|\hat f\|^2_{L^1(\mathbb{R}^d)}}/{(KL)}$ of the generalization error for random Fourier features with one hidden layer and the same total number of nodes $KL$, in the case the $L^\infty$-norm of $f$ is much less than the $L^1$-norm of its Fourier transform $\hat f$. This understanding of an optimal distribution for random features is used to construct a new training method for a deep residual network. Promising performance of the proposed new algorithm is demonstrated in computational experiments.

In quantitative finance, modeling the volatility structure of underlying assets is vital to pricing options. Rough stochastic volatility models, such as the rough Bergomi model [Bayer, Friz, Gatheral, Quantitative Finance 16(6), 887-904, 2016], seek to fit observed market data based on the observation that the log-realized variance behaves like a fractional Brownian motion with small Hurst parameter, $H < 1/2$, over reasonable timescales. Both time series of asset prices and option-derived price data indicate that $H$ often takes values close to $0.1$ or less, i.e., rougher than Brownian motion. This change improves the fit to both option prices and time series of underlying asset prices while maintaining parsimoniousness. However, the non-Markovian nature of the driving fractional Brownian motion in rough volatility models poses severe challenges for theoretical and numerical analyses and for computational practice. While the explicit Euler method is known to converge to the solution of the rough Bergomi and similar models, its strong rate of convergence is only $H$. We prove rate $H + 1/2$ for the weak convergence of the Euler method for the rough Stein-Stein model, which treats the volatility as a linear function of the driving fractional Brownian motion, and, surprisingly, we prove rate one for the case of quadratic payoff functions. Our proof uses Talay-Tubaro expansions and an affine Markovian representation of the underlying and is further supported by numerical experiments. These convergence results provide a first step toward deriving weak rates for the rough Bergomi model, which treats the volatility as a nonlinear function of the driving fractional Brownian motion., Comment: 28 pages, 4 figures

We provide an overview of the methods that can be used for prediction under uncertainty and data fitting of dynamical systems, and of the fundamental challenges that arise in this context. The focus is on SIR-like models, that are being commonly used when attempting to predict the trend of the COVID-19 pandemic. In particular, we raise a warning flag about identifiability of the parameters of SIR-like models; often, it might be hard to infer the correct values of the parameters from data, even for very simple models, making it non-trivial to use these models for meaningful predictions. Most of the points that we touch upon are actually generally valid for inverse problems in more general setups.

We develop a data-driven methodology based on parametric It\^{o}'s Stochastic Differential Equations (SDEs) to capture the real asymmetric dynamics of forecast errors. Our SDE framework features time-derivative tracking of the forecast, time-varying mean-reversion parameter, and an improved state-dependent diffusion term. Proofs of the existence, strong uniqueness, and boundedness of the SDE solutions are shown under a principled condition for the time-varying mean-reversion parameter. Inference based on approximate likelihood, constructed through the moment-matching technique both in the original forecast error space and in the Lamperti space, is performed through numerical optimization procedures. We propose another contribution based on the fixed-point likelihood optimization approach in the Lamperti space. All the procedures are agnostic of the forecasting technology, and they enable comparisons between different forecast providers. We apply our SDE framework to model historical Uruguayan normalized wind power production and forecast data between April and December 2019. Sharp empirical confidence bands of future wind power production are obtained for the best-selected model., Comment: 28 pages and 11 figures

In this paper, we consider the filtering problem for partially observed diffusions, which are regularly observed at discrete times. We are concerned with the case when one must resort to time-discretization of the diffusion process if the transition density is not available in an appropriate form. In such cases, one must resort to advanced numerical algorithms such as particle filters to consistently estimate the filter. It is also well known that the particle filter can be enhanced by considering hierarchies of discretizations and the multilevel Monte Carlo (MLMC) method, in the sense of reducing the computational effort to achieve a given mean square error (MSE). A variety of multilevel particle filters (MLPF) have been suggested in the literature, e.g., in Jasra et al., SIAM J, Numer. Anal., 55, 3068--3096. Here we introduce a new alternative that involves a resampling step based on the optimal Wasserstein coupling. We prove a central limit theorem (CLT) for the new method. On considering the asymptotic variance, we establish that in some scenarios, there is a reduction, relative to the approach in the aforementioned paper by Jasra et al., in computational effort to achieve a given MSE. These findings are confirmed in numerical examples. We also consider filtering diffusions with unstable dynamics; we empirically show that in such cases a change of measure technique seems to be required to maintain our findings., Comment: 48 pages