Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
66,531 characters · 14 sections · 96 citation commands
Fast Estimation of Bayesian State Space Models Using Amortized Simulation-Based Inference The views expressed in the paper are solely those of the authors and do not necessarily represent the official position of the Bank of Russia. The Bank of Russia is not responsible for the contents of the paper.,We are grateful to Dmitry Gornostaev, Sergey Ivaschenko, Petr Milyutin, Denis Shibitov and participants of the 15th International Conference on Computational and Financial Econometrics (CEF 2021), XXIII Yasin (April) International Academic Conference on Economic and Social Development and the 2nd International Conference on Econometrics and Business Analytics (iCEBA) for their helpful comments and suggestions.,
Bayesian state space models are widely used in applied macroeconomics. They are so widespread due to the fact that many macroeconomic and econometric models can be written in the form of state space models for subsequent estimation on real data. For example, various kinds of filters and semi-structural filters (see hodrick1997postwar, laubach2003), models with stochastic volatility (see kim1998stochastic, justiniano2008time, carriero2019large), time-varying models (see 10.2307/1912559, primiceri2005time, korobilis), mixed frequency models (see chiu2012estimating, schorfheide2015, schorfheide2021real), dynamic factor models (see otrok1998bayesian, stock2011dynamic), dynamic stochastic general equilibrium models (see smets2003estimated,smets2007shocks, FERNANDEZVILLAVERDE2016527) and agent-based models (see lux2018estimation, gatti2020rising). Bayesian parameter estimation makes it possible to mitigate the lack of long time series \oldFootnote{In this paper, we focus on time series models, but the proposed algorithm can easily be transferred to hidden space models of other types with minor modifications in the architecture.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi.
Despite their flexibility, in practice, estimating Bayesian state space models is a rather difficult task. Sampling algorithms are based on an iterative sampling scheme for model parameters and states and rely on steps such as Gibbs sampling (see casella1992explaining), Metropolis-Hastings (see chib1995understanding), Hamiltonian Monte Carlo (see neal2012bayesian) or sequential Monte Carlo (see del2006sequential). Even in cases where the model is linear and Gaussian or discrete with respect to states, sampling can take from tens of minutes to hours. In systems of a more general form, researchers use particle filters (see andrieu2010particle, chopin2013smc2), as a result of which estimation time only increases. Sampling algorithms are exact in the sense that they converge to the posterior distribution as the number of iterations tends to infinity, although the number of iterations required for an estimate close to the posterior distribution can be large. An alternative to them are optimization algorithms. The most common of them in the context of state space models is the variational Bayes algorithm \oldFootnote{See Chapters 3 and 5 of beal2003variational and gunawan2021variational for examples of using the Bayesian variational algorithm in estimating of state space models.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi (see wainwright, hoffman2013stochastic). The variational Bayes algorithm relies on minimizing the Kullback-Leibler (or any other) divergence between the approximation and the true posterior distribution. It is often faster than sampling algorithms, but its running time is also large, especially in cases where the optimization steps cannot be written in analytical form.
In this paper, we propose a fast algorithm for estimating Bayesian state space models that is based on the principles of simulation-based inference (see cranmer2020frontier). The speed of the algorithm is achieved by amortizing the task of constructing the posterior distribution of states, that is, by pretraining a model (a neural network, in our case) that predicts its posterior distribution from the data. In doing so, we focus only on states and do this for two reasons. First, in many problems it is the states, not the model parameters, that are of particular interest. For example, in the detrending problem (see orphanides2002unreliability), the trend and cycle components, which are determined by hidden states, are of main interest. Second, simulation-based inference parameter estimation has been investigated in many other papers (see Appendix A of lueckmann2021benchmarking) and can be easily combined with the approach discussed in this paper. The problem of constructing posterior distribution of hidden states is much more complicated due to its dimensionality and has not been studied much in the literature.
Section (ref) describes the algorithm for estimating the posterior distribution of the model. Section (ref) is devoted to the study of application, practical characteristics and comparison of the algorithm with commonly used alternatives. Section (ref) describes related work. Section (ref) discusses issues that have not been included in the paper but are important in the context of the proposed algorithm. The conclusion is presented in Section (ref).
Our methodology for estimating state space models is based on the idea of estimating Bayesian parameters proposed by beaumont2002approximate, blum2010non and developed by papamakarios2016fast. The essence of the methodology is to simulate the joint distribution of parameters and data, and then predict the parameters conditional on data.
Formally, during the first step, a dataset is simulated from a model with a prior distribution $p(\theta)$ and a likelihood function $p(y | \theta)$. The $i$-th point consists of parameters $\theta_i \sim p(\theta)$ and observed variables $y_i \sim p(y | \theta_i)$. In the second step, an estimator is fitted to predict the distribution $p(\theta | y)$ or its characteristics, for example, by minimizing the cross-entropy between the simulated data and some parametric family of distributions $q_{\varphi} (\theta|y) p(y)$:
where $\varphi$ is the vector of parameters of distribution $q_{\varphi} (\theta | y)$, $N$ is the number of simulations. We will omit $log\,p(y_i)$ term due to its independence from $\varphi$ in the following.
The estimated distribution will tend to the posterior for any $y$ with an infinite number of simulations and a sufficiently flexible parametric family $q_{\varphi}$ (see Appendix (ref) for an informal proof). This means that once estimated, the conditional distribution $q_{\varphi^{*}} (\theta | y) $ can be used for any data, does not require re-estimation of the model and can be calculated almost instantly.
As can be seen, this method for estimating posterior distributions (hereinafter, we will call it NPE following papamakarios2016fast) does not require knowledge of $p(\theta)$ and $p(y|\theta)$ in an explicit form but relies only on the ability of simulating data, which is natural for most models. So, it falls into the category of simulation-based inference (SBI) methods (or likelihood-free inference).
The parameters of state space models can be estimated using the procedure presented in Subsection (ref), however the goal of this paper is to estimate the hidden states. It is easy to see that if we replace the parameters with hidden states in the loss function, then the procedure described above remains valid (see Appendix (ref)). Thus, in general, the algorithm for finding the posterior distribution can be written as shown below:
\,\,\,\,\,\, Despite the simplicity of Algorithm (ref), there are several practical difficulties when moving from finding posterior distribution for the parameters to finding the distribution of states. First, it is a large dimension of the hidden space. Despite the presence of various kinds of flow transformations (see rezende2015variational), which are often used in SBI and in many cases recover joint distributions adequately, their application for problems of this size, complicated by amortization, is costly to compute and is associated with optimization challenges. Second, it is the large dimensionality of the data. It is nearly impossible to find summary statistics that reduce the dimensionality of the data for hidden states, in contrast to the problem of parameter estimation, where this is common practice (see, for example, SIR (T.9) and Lotka-Voltera (T.10) models in lueckmann2021benchmarking). To overcome these problems and simplify the task of training the model, we avoid modeling the dependencies between variables focusing on characteristics of marginal distributions and introduce an inductive bias for $q_{\varphi} (s| y)$.
Moving from the joint to the marginal distributions is equivalent to subdividing the task into a set of one-dimensional tasks. In such case, the log-probability for the parametric family $q_{\varphi} (s | y)$ is written as:
where $t$ and $k$ represent the time period and the state index in state vector. Note that for each state, the vector of parameters $\varphi_{t,k}$ is generally its own and $\varphi$ consists of a set of these vectors.
We use the normal distribution with mean $m_{\varphi_{t,k}} (y)$ and standard deviation $\sigma_{\varphi_{t,k}} (y)$ for $q_{\varphi_{t,k}} (s^{t,k} | y)$. It is easy to show that for sufficiently flexible $m_{\varphi_{t,k} } (y)$ and $\sigma_{\varphi_{t,k}} (y)$ such an approximation exactly recovers the mean and standard deviation of the true posterior distribution (see Appendix (ref)).
Moreover, the choice of normal distribution and cross-entropy as a loss function can be relaxed. Therefore, a mixture of normals or a small-scale flow-based model can replace normal distribution for marginal densities and any M-estimator (see Chapter 5 in van2000asymptotic), such as quantile regression loss, can be used instead of cross-entropy.
A natural architecture that is similar to filtering and smoothing in state space models is a Bidirectional Recurrent Neural Network (Bidirectional RNN). By sharing parameters for models with a «near» stationary data generation process \oldFootnote{We use the term «near» stationary data generation process to emphasize the possibility of using different state distributions for the initial time period.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi \oldFootnote{One can always connect several networks, if there is a shift in the data generation process.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi, this structure allows to move away from estimating the neural network for each state separately and calculate the loss function for all states in a single pass over the data.
To make the network architecture more flexible, we use data convolution of various lengths as RNN input in addition to the raw data, and transform the output of RNN by applying a linear transformation or a fully connected neural network.
To illustrate the properties of the proposed method, we estimate three models: stochastic volatility model, non-linear DSGE model and seasonal adjustment model with structural breaks in seasonality.
The stochastic volatility model (see kim1998stochastic) is a classic example for testing various states estimation algorithms in state space models (see tan2020conditionally). The model specification is exactly the same as tan2020conditionally and presented in Appendix (ref). To simplify the learning process, the logarithm of absolute values of the observed data is fed into the neural network as an input (details of the network architecture and the training algorithm are presented in Appendix (ref)). The model is trained on 20,000,000 generated series with a length of 800 to 1,200 observations.
First, we demonstrate the quality of the algorithm on the data sampled from the data generation process. Figure (ref) shows examples of the true values of volatility logarithms and mean of posterior distribution approximation ($\pm$ 2 standard deviations). The figure demonstrates that the true values are quite well estimated by the neural network. As a benchmark for future studies, we also present the negative log-likelihood (hereinafter, NLL) and mean squared error (hereinafter, MSE) on a randomly generated 1,000,000 runs in Table (ref).
Unfortunately, calculation of the accuracy metrics for other algorithms, such as MCMC or stochastic variational Bayes, is difficult to compute (see the discussion on quality metrics for SBI in lueckmann2021benchmarking). Our focus therefore lies on the comparison of the results for NYSE and GBPUSD datasets similarly to tan2020conditionally \oldFootnote{Calculating metrics, such as C2ST, are not informative for the joint distribution. Information about correlations is important for classification, but the algorithm used assumes the diagonal covariance matrix. Calculation of C2ST for marginal distributions requires training the number of classifiers equal to the number of hidden states. So, we use only visual analysis.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi. The results are compared with the adaptive MCMC algorithm based on the mixture of normals for chi-square distribution approximation, which was proposed by kim1998stochastic, and the stochastic variational Gaussian approximation (hereinafter, VB) with a sparse precision matrix (both algorithms are given in Appendix (ref)). As can be seen in Figure (ref), although the neural network estimates are a bit noisy, they are close to the MCMC algorithm, which serves as the gold standard, as well as the variational Bayes algorithm, which is one of the fastest and most accurate approximations. It should be noted that the NYSE dataset is almost twice the maximum size of the simulation, and the neural network has never seen data of that length. Nevertheless, the trained model copes with this task.
The stochastic volatility model, while contains many hidden states, is a univariate model, both in terms of data and output \oldFootnote{The latter could potentially be an advantage, though, as information from different sources can help in training neural network parameters. Moreover, it is possible to train a univariate model for each dimension of state vector, which reduces the problem to the previous one in terms of output.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi. The DSGE model was chosen to test the amortized SBI algorithm in multivariate context. This class of models is widely used by macroeconomists, both for practical purposes (see linde2016challenges) and in academic research (see walsh2010). Although solving non-linear DSGE models is beyond the scope of this paper, we seek to demonstrate that the proposed algorithm works well for models where filtering and likelihood estimation cannot be executed using the Kalman filter as in the case of linear models \oldFootnote{Log-linearized versions of these models are often used to overcome the computational difficulties with solving and estimating non-linear models.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi. A simplified DSGE model \oldFootnote{We exclude the inflation target shock from the model and inflation expectations from the observed variables to avoid the issues of missing variables and their impact on the result. A number of additional experiments have shown that using trained values for neural network inputs in place of missing variables (see lueckmann2017flexible) and introducing additional dummy variables to the RNN input can cope with this task. However, we will focus here on a simpler version of the model to separate the effect of multiple observed variables from the effect of missing data. Estimates of DSGE models using SBI will be the subject of a separate paper where we will also touch on this issue.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi with stochastic volatility from diebold2017real was chosen for these reasons. This model can be solved using standard algorithms for linearized models (see blanchard1980solution, anderson1985linear, klein2000using, sims2002solving). Nonlinearity is introduced after the solution step as the time-varying volatility of model shocks. The neural network is estimated on 50,000,000 generated datasets with a length of 180 to 200 points and compared with the adaptive MCMC algorithm (model description, neural network architecture and MCMC implementation are described in (ref)).
To illustrate the properties of NPE, we concentrated on the estimation of stochastic volatilities and present graphs and metrics for these states. However, unobservable shocks (more precisely, the logarithms of their absolute values) were also used in the estimation as a hint on the intermediate outputs of the neural network. Figure (ref) shows that the trained neural network results are similar to MCMC for the US data from 1964Q2 to 2011Q1 (latest vintage in diebold2017real). As for the previous model, Table (ref) shows the NLL and MSE on 1,000,000 randomly generated datasets for future comparisons.
In addition to problems based on well-verified formulas for transition and observation equations, SBI is also suitable for those models where simulations are the primary focus. Usually, this occurs in models where deriving equations is too cumbersome or simply impossible to compute. Special cases are tasks where it is easy to generate a lot of different data which shows the model how it should behave in various situations.
As an example, we show how a quarterly seasonal adjustment model that considers structural shifts in the seasonal component can be built. As will be shown below, a traditional X13 ARIMA-SEATS procedure (see x13arimaseats) does a poor job of this. An additional advantage is the automatically generated credible intervals.
(ref) presents the procedure for generating artificial series with a length of 40 to 80 quarters. In fact, it consists of generating a seasonal and non-seasonal component with the probability of a shift appearing in the seasonal part. Thus, the resulting series may not contain a break.
We compare NPE not with a sampling algorithm, but with X13 for measuring the quality, in contrast to the previous two models. The purpose of this experiment is to demonstrate how one can easily generate examples of model behavior, thus specifying an implicit Bayesian model. In practice, it is usually difficult to construct a fast MCMC algorithm in such cases. However, comparison with other algorithms that solve the same practical problem is of interest from the point of view of estimating the performance of the proposed algorithm.
Figure (ref) shows random examples illustrating how the proposed procedure and X13 behave on series with a shift in seasonality (the gray area shows 1.5 years around the shift). X13 does not adequately cope with the task of detecting seasonality around the quarter of shift, while NPE, on the contrary, is robust. We calculated the MSE on 100,000 randomly generated runs for X13 and on 1,000,000 runs for the neural network. The error for NPE is smaller, as can be seen from Table (ref). The smaller errors in comparison to X13 are not unexpected by themselves since cross-entropy optimization should provide the estimator with the lowest MSE on this dataset. However, the gap between the errors on the series with and without shifts further emphasizes that the algorithm proposed in the paper is more accurate than the alternative widely used among macroeconomists.
As noted above, NPE works almost instantly due to possessing an amortization property. Depending on the task, calculating an approximation of the posterior distribution on the CPU (Intel(R) Core(TM) i7-8750H CPU @ 2.20GHz, 16GB RAM) takes tenths of a second for a pretrained model. Estimation of neural network parameters on Pytorch \oldFootnote{Generation of artificial data for the model is run on the CPU and training on the GPU.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi (see NEURIPS2019_9015) using GPU (NVIDIA GeForce RTX 2070) takes about 12, 12 and 2 hours for the stochastic volatility model, DSGE model and seasonal adjustment model, respectively. Comparing the running time for amortized algorithms with alternatives that do not possess such properties comes with difficulties. On the one hand, Bayesian model estimation for a fixed dataset takes less time in our examples when MCMC or VB algorithms are used (see Table (ref)) \oldFootnote{We choose the number of iterations based on convergence of parameters for MCMC (half of iterations is a burn-in period) and convergence of loss for VB.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi. This is because the amortized algorithm requires pretraining of a neural network. On the other hand, a neural network trained once can be used for various datasets (including data of different lengths) which is an advantage in the case of multiple estimations.
The implementation of the NPE algorithm uses already available neural network libraries (coding a neural network architecture usually takes only tens of minutes since the layers of the neural network are already implemented in the respective libraries) and almost always turns out to be much easier than the implementation of MCMC and even sometimes variational estimation. The implementation of MCMC requires the derivation of the sampling algorithm and the writing of program code, which is often much more complicated than for NPE. Like NPE, stochastic VB algorithms are just as easy to implement in most cases (if automatic differentiation packages are used and there are no modules that require differentiation coding). The key difference is that the joint density of parameters, hidden states and data is used instead of data sampling procedure.
An important point for stochastic optimization algorithms (NPE and VB) is also the ability to almost effortlessly transfer the computation to the GPU.
The algorithm proposed in this paper is closely related to several directions presented in the literature. Our work is part of the literature on SBI (or likelihood-free inference) algorithms (see cranmer2020frontier). Until recently, Bayesian direction in this field has developed mainly as approximate Bayesian computations (hereinafter, ABC). ABC approximate the likelihood functions by introducing an auxiliary likelihood function that depends on the distance between summary statistics of the data and the corresponding statistics of the simulated data. Classic sampling algorithms are then used (see sisson2018handbook). The progress of machine learning algorithms, and in particular neural networks, has given rise to a whole family of algorithms that directly train posterior distributions (see papamakarios2016fast, lueckmann2017flexible, greenberg2019automatic, durkan2020contrastive), likelihood functions (see wood_statistical_2010, lueckmann2018likelihood, brehmer2020mining, papamakarios2019sequential) or likelihood function ratios (see brehmer2020mining, hermans2020likelihood, durkan2020contrastive). These methods are not sensitive to tolerance hyperparameter and chosen distance between simulation data and summary statistics, contrary to ABC. However, to the best of our knowledge, apart from a few papers on probabilistic programming (see le2017inference, baydin2019etalumis, munk2019deep), researchers concentrate mainly on parameters, not states. Research on probabilistic programming has two key differences from the approach proposed here. First, it uses the pretrained neural network as a proposal distribution for importance sampling, rather than directly to approximate the posterior distribution. Second, a neural network in probabilistic programming takes into account the relationship between variables to achieve smaller variance for the importance sampling weights, which is a considerably more difficult task in terms of optimization \oldFootnote{It is also worth noting that probabilistic programming uses state sampling which depends on the sampled states of the previous period. This can lead to accumulation of approximation errors over long periods. Such an architecture is not quite suitable for direct approximation. However, this is not critical for subsequent resampling, especially if sequential importance sampling is used instead of importance sampling.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi. Moreover, unlike this paper, implementation from scratch or modification of probabilistic programming algorithms for tasks that do not fit into the framework of standard libraries \oldFootnote{See, for instance, PyProb.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi is quite complicated, since it requires a deep knowledge of the addressing of random variables.
Our research is also closely related to the estimation of economic models through simulations. The simulated method of moments and its modifications (see mcfadden1989method, 10.2307/2951768, gallant1996moments) are common in the frequentist estimation of the structural parameters of models \oldFootnote{See a list of applications in carrasco2002simulation.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi. There is a similar field of research that estimates model parameters based on minimizing various divergences between simulated and real data (see nickl2010efficient, kaji2020adversarial). gallant2009determination \oldFootnote{ gallant2013generalized also developed a version of the Bayesian simulation method of moments based on a particle filter for models with hidden states.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi proposed a Bayesian version of the simulated method of moments. A recent paper by fen2022fast uses sequential (non-amortized) NPE for Bayesian parameter estimation. As for SBI, not many papers devoted to simulated estimation of states rather than parameters exist. The closest known to us is the paper by gatti2020rising, where the authors estimate states (output and investment gap forecasts) on artificial data using nonparametric kernel estimation. The idea is very close to SBI and to what is proposed here, but it is computationally difficult with a large number of hidden states as the authors themselves note.
Meta-learning (see finn2019meta) is close to SBI in its mathematical formulation. Like SBI, meta-learning is based on the idea of learning from many similar tasks (see vinyals2016matching). The key differences are purpose and data. Unlike SBI, meta-learning focuses on the task of predicting rather than finding the posterior distribution. Furthermore, meta-learning usually works with real data, not simulated ones.
Many works on variational autoencoders have been devoted to the amortization of finding the distribution of hidden states (see kingma2019introduction). However, there are a number of differences from this paper. First, the data generation process is usually specified with a rather flexible model such as a neural network (see kingma2013auto) or a Gaussian process (see dai2015variational) rather than more classical models where the hidden states have greater identifiability and interpretability. In addition, the model is usually non-Bayesian in nature \oldFootnote{Basically, neural networks are used as a data generation model, which by their nature are frequentist. Moreover, despite the fact that many applications use dropout for regularization (srivastava2014dropout), which has a Bayesian interpretation (see kingma2015variational and gal2016dropout), the parameters are common to all data. This is ideologically different from the idea of amortizing models.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi. Second, the loss function that is minimized is the KL divergence between the approximate posterior and posterior distributions, while in SBI it is the KL divergence between the posterior and approximate posterior distributions. The asymmetry of KL divergence leads to the fact that, with few exceptions (see tran2017hierarchical), there are not enough simulations to train variational autoencoders and one must calculate the probabilities of the data generation process. Also, in the case of diagonal approximation (as in Section (ref)), this leads to underestimation of variance (see blei2017variational). Third, real, not artificial, data are used for training variational autoencoders as for meta-learning.
As has been shown in many papers (see lueckmann2021benchmarking), amortization leads to the need for longer training of SBI algorithms than their sequential counterparts. Despite this, we use the amortized NPE algorithm for two reasons. First, sequential SBI algorithms are usually applied to the problems of small dimension (the output of neural network dimension), and their adaptation to the high-dimensional problem of finding the posterior distribution of states is not trivial and requires solving more practical issues. In particular, if one tries to focus on marginal densities, as is done in this paper, the states generated for new rounds will not look like posterior distribution due to the lack of dependence between variables. The states will be quite noisy, which worsen convergence in most cases. Secondly, in contrast to some research on SBI, we do not set the task to finding the best algorithm under the constrained budget for the number of simulations (see lueckmann2021benchmarking). Sequential algorithms give a significant gain for such tasks. However, the main goal of this paper is to build an algorithm that replaces long-running alternatives (as in the first two examples) or helps to estimate models where other algorithms fail (as in the third example) \oldFootnote{It is implicitly assumed that a large number of model simulations can be performed in adequate time.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi . Amortization is a great property that helps to solve this problem if the model is frequently re-estimated.
Many issues related to the estimation of the posterior distribution of states are beyond the scope of this paper and require further research. Some of them are discussed below.
The posterior distributions estimated using the proposed algorithm, although close to the MCMC results, nevertheless differ slightly. The results are slightly noisy when estimating stochastic volatility, while in the DSGE model they are biased. This signals an opportunity for further improvements of the neural networks by increasing the flexibility of the neural network architecture, the number of simulations or by modifying the training procedure. It is well known that for a certain learning rate schedule, stochastic optimization procedures converge to one of the local optima in the asymptotics (see Chapter 5 in kushner2003stochastic). The exact (even local) optimum is not achieved with a finite number of iterations \oldFootnote{This usually means that the learning rate does not tend to zero.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi \oldFootnote{ mandt2017stochastic provide intuitions about the behavior of the estimation procedure at non-zero learning rates.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi. An insufficiently flexible network and/or a small number of observations in the neighborhood of real data can lead to a situation where the model is unable to predict the posterior distribution accurately, even at the optimum \oldFootnote{A good example of such type an improvement in the field of text analysis is the GPT-3 model (see brown). It has reached a fundamentally new level compared to previous models due to an order of magnitude more parameters than previously used and a huge dataset.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi. Moreover, the quality of the posterior distribution approximation is likely to deteriorate with increasing problem dimensionality. Hence, one of the main tasks for the future is to study the relationship between scalability, approximation quality and neural network training time.
The mean-field Gaussian approximation considered here is usually not a problem from a practical point of view, because in most cases, researchers are interested in the first and second moments of the marginal distributions of states. Although extensions have clear theoretical solutions (M-estimators for estimating other characteristics and more flexible families of distributions), their practical implementation requires further research \oldFootnote{In a number of preliminary experiments that were not included in the paper, we saw that the quantile loss also shows good results for the marginal distributions.}\futurelet\nextToken \ifx\footnote\nextToken\textsuperscript{,}\fi.
We have bypassed the issues of forecasting and missing variables, which are related in the sense that the forecasting problem can be thought of as a problem of constructing a posterior distribution for the missing variables on the forecasting horizon. To deal with missing variables, models can be extended by introducing additional dummy variables as one of the inputs of the neural network, showing the presence of a miss, and/or by filling in the miss with learnable parameters as done in lueckmann2017flexible. A similar method or alternatives based on meta-learning ideas (see harrison2018meta), where only the predicted variables are used as the neural network output, can be applied to build prediction models.
NPE has both advantages and disadvantages in terms of speed as shown in Section (ref). Therefore, the choice to use NPE or not should depend on the situation. We recommend using the NPE algorithm if frequent re-estimation of the model is expected, or if alternative algorithms are slow, or fail to do the job at all. At the same time, we also advise to verify the trained algorithm before use by comparing it with alternative ones for approximating the posterior distribution (when they are not too slow). If such verification is impossible, it is recommended to carry out at least a visual analysis on artificially generated data.
The amortized simulation-based algorithm proposed in this paper for estimating hidden states of Bayesian state space models provides an alternative to already existing algorithms in this field. In contrast to many previous papers, we consider a new approach that approximates posterior marginal distribution of states and that does not rely on probability density functions for prior distributions, transition and observation equations, but that uses only simulations of artificial data.
The NPE algorithm shows results similar to other algorithms for the stochastic volatility and DSGE models but after training, it works nearly instantly. In addition, as shown in the example with seasonal adjustment, it also performs well on tasks where the Bayesian model is not specified directly but rather through the process of simulating various situations and correct behavior in them.