{'fraction_non_alphanumeric': 0.03933860270736804, 'fraction_numerical': 0.009630235304806032, 'mean_word_length': 4.393189612934836, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the Simulation of Polynomial NARMAX Models

In this paper, we show that the common approach for simulation non-linear stochastic models, commonly used in system identification, via setting the noise contributions to zero results in a biased response. We also demonstrate that to achieve unbiased simulation of finite order NARMAX models, in general, we require infinite order simulation models. The main contributions of the paper are two-fold. Firstly, an alternate representation of polynomial NARMAX models, based on Hermite polynomials, is proposed. The proposed representation provides a convenient way to translate a polynomial NARMAX model to a corresponding simulation model by simply setting certain terms to zero. This translation is exact when the simulation model can be written as an NFIR model. Secondly, a parameterized approximation method is proposed to curtail infinite order simulation models to a finite order. The proposed approximation can be viewed as a trade-off between the conventional approach of setting noise contributions to zero and the approach of incorporating the bias introduced by higher-order moments of the noise distribution. Simulation studies are provided to illustrate the utility of the proposed representation and approximation method.


I. INTRODUCTION
In the field of data-driven modelling, prediction error has been the primary choice of identification criterion for many decades. This choice is well-justified due to its ease of use and the maximum-likelihood interpretation of the one-step ahead predictor with respect to the sum-of-squared error identification criterion [1]. This interpretation holds true under the assumption that the underlying noise process is zero-mean and that the true data-generating system can indeed be represented by the chosen model class. It has been observed (for example, in [2], [3]) that the presence of an auto-regressive component in the prediction model can allow a relatively simpler model to achieve optimistic prediction results for an otherwise complex dynamical system. This observation leads to the following question -given a set of data measured from the unknown dynamical system and an identified model, how can one ensure that the model not only predicts the measured data, but also captures the dynamical structure of the true system. This question has been long recognized in literature (c.f. [4]). A key challenge that is inherent to this question is that one should be able to assess the quality of the identified model across multiple model classes, for example, to compare the quality of a linear model with that of a non-linear model.
Model validation is an important ingredient of the system identification framework. A common approach in model validation is to order models based on a pre-defined measure of complexity, and then add levels of complexity in a model until the data can no longer invalidate the model [5]. This approach is in line with the Occam's Razor principle that suggests that one should choose the simplest model that explains the data. However, the approach becomes quite cumbersome and computationally infeasible when identification is done in a black-box setting that spans across multiple model classes. Another common approach to model invalidation is to estimate the generalization error of an identified model (i.e. the performance of an identified model on an independent data-set). When a sufficiently large amount of data is available, the user can set aside a part of the data that is not used during model estimation. The independent dataset can then be used to estimate the generalization error of all models identified based on the data used for estimation. This method is called cross-validation [6]. When the amount of data is insufficient, information criterion such as AIC (Akaike Information Criterion) or BIC (Bayesian Information Criterion) can be used to estimate the generalization error.
An alternative approach that deals with this challenge is the control of complexity by using regularization. In these methods, the cost function is modified to include a penalty term that controls the complexity of the identified models as part of the optimization problem. However, in a black box setting that spans multiple model classes, finding a suitable regularization penalty is challenging.
A third alternative, and one that is of interest in this contribution, is to make use of simulation error to judge the quality of an identified model. Simulation models do not make use of the past output measurements, as opposed to prediction models. As a result, simulation error is typically more sensitive to mismatch between the model and the system structure than prediction error (see [2], [7]). Moreover, simulation models can be defined for models belonging to different model classes. Simulation error has been used in the past for non-linear structure selection, for example in [2].
The notion of a simulation model can be understood in multiple ways. A simulation model can be viewed as an infinite horizon prediction model. This effectively negates the auto-regressive components of a prediction model. A simulation model can also be interpreted as the "deterministic response" of the system. As a result of this interpretation, simulation models are often reported in the literature as prediction models with noise terms set to 0 [8, ch. 5]. In the case of linear systems, it is always possible to lump the noise sources as an additive term on the output of the deterministic response of the system. Hence, the approach of setting noise contributions to 0 can be intuitively linked to the idea of obtaining the "deterministic response" of the system. However, in the non-linear setting, neglecting the process noise terms leads to a biased response (see [9], [10]).
In this paper, two contributions are made with respect to the computation of simulation models for polynomial NARMAX models. In Section IV-A, a Hermite polynomialbased representation is proposed. The proposed representation provides a convenient way to link prediction models and simulation models, under the assumption that the noise is distributed normally. In Section IV-B, it is shown that computing a simulation model for a finite order NARMAX model may lead to an infinite order NIIR model. For such cases, an approximation method is proposed. Numerical illustrations are provided in V.

II. PROBLEM DESCRIPTION
Consider the discrete-time single-input single-output (SISO) polynomial NARMAX model form: where, f (.) is a polynomial non-linear function, u k , y k ∈ R are the input and the output of the system at time instant k, ξ is a zero-mean white noise process assumed to be independent of u, and n u , n y and n ξ are the maximum time lags for the respective signals. In the sequel it is assumed that f in (1) is Bounded Input Bounded Output (BIBO) stable with respect to the deterministic input u and that the initial conditions {u k , y k , ξ k } for k ≤ 0 are 0. Note that (1) is not the most general representation of non-linear systems, nonetheless, a large family of systems can be modelled within this class.
Given a model of form (1), one possible approach to compute the corresponding simulation model is to set the noise contributions to 0. This approach works well when the function f is linear because it is always possible to separate the noise contributions in the function f to form an additive noise model of the form where g and h are linear functions. If a simulation model is described as the deterministic response of the system, then it can be expressed as the conditional expectation of the output E[y k ] taken with respect to the distribution of the noise. Since it is assumed that the noise process is zero-mean, the simulation response y s can be written as y s,k := E[y k ] = g(y s,k−1 , . . . , y s,k−ny , u k , . . . , u k−nu ).
(3) Notice that this is the same as setting all noise contributions in (1) to 0. However, when f is a non-linear function, it is not trivial to separate the noise contributions into a separate noise process that appears additively in the output of the system. Furthermore, in order to compute the conditional expectation of (1), one must take the expectation over terms that include polynomial factors of the past outputs y and noise ξ. In order to compute such a conditional expectation, one must take in to account the higher order moments of the noise process and the non-linear dependencies of the output terms y on the noise terms ξ. As a result, the simple approach of setting the noise terms to 0 results in a biased simulation response.
Continuing the line of reasoning from the linear case, one can represent the simulation response y s of (1) as: where the conditional expectation is taken with respect to the probability density p ξ of ξ and the input sequence The problem can be stated as follows. Given an input-output sequence {u k ; k ∈ [1, N ]}, and a model of form (1), compute the simulation model, as given in (4).
It is important to note that in order to compute simulation models as in (4), higher order moments of the noise distribution must be known. Hence, in the sequel, we assume that the noise process ξ k is i.i.d and is distributed as and hence, can be included in the model formulation in (1).

III. MATHEMATICAL BACKGROUND
Under the given assumptions that the noise process is distributed normally, the noise process terms ξ k that appear in model (1) may be transformed to Hermite polynomials [11,Ch. 18] in ξ k . For completeness, a brief overview of the relevant properties of Hermite polynomials is presented.
Let (R, B R , p x ) be a probability space on the real line R, where B R is the Borel σ-algebra on R, and p x is the probability measure of the Gaussian distribution. Consider the Hilbert space L 2 (R, B R , p x ) of random functions f that satisfy, where E is associated with the inner product for any f, g ∈ L 2 (R, B R , p x ). Definition 1: Hermite polynomials He n (x) of degree n are defined as Hermite polynomials form a closed and complete orthogonal basis in the L 2 space with the inner product defined in (7) (see [11]), i.e., where δ n,m = 1 when n = m and 0 otherwise. Hermite polynomials possess a number of properties that will be used in order to derive simulation models automatically from a given stochastic model of form (1). Since Hermite polynomials form a complete orthogonal basis in L 2 (R, B R , p x ) with inner product defined in (7), any function f ∈ L 2 (R, B R , p x ) can be expressed as a convergent linear combination of basis functions Given a monomial in x ∈ R, one can compute an equivalent representation in terms of Hermite polynomials He(x) using the following inverse relationship where ⌊.⌋ is defined as the floor operator, i.e., ⌊n⌋ = q, where q ∈ Z is the largest integer such that q ≤ n. Hermite polynomials of normally distributed random variables have useful properties. Let X be a random variable distributed as N (0, 1). It can be shown that The following result can be obtained by using (9) It has been shown in [12] that for X 1 , X 2 ∼ N (0, 1), In particular, if X 1 , X 2 are independent, we get
It should be noted that an equation of form (16) can be easily transformed to the equivalent simulation model as in (18) without the use of Hermite polynomials, since the moments of a Gaussian distribution are known. However, the use of Hermite polynomials yields an equivalent prediction model, as in (17), that offers a convenient representation. It turns out that, in order to obtain a simulation model like (18), one must only "switch off" a number of Hermite polynomial terms in (17), as per Equations (12), (13) and (15)). An example will be used to illustrate the method to compute the transformation, followed by a general result.
Example 1: Consider a system described by the following equation where {u k } and {ξ k } are i.i.d sequences distibuted as N (0, 1). The given equation can be re-written as follows Now, the simulation model can be computed by taking the expectation of (20) with respect to the noise process ξ, Observe that if the simulation model was computed by setting the contributions of the noise terms to 0, one would obtain the following simulation model The two simulation models are offset by a scalar factor 1.
In representation (17) (and (20)), one must only set all nonzero order Hermite polynomials of the noise term to 0 in order to obtain a simulation model as in (18) and (21). For a general result, the following equivalent representation of (16) will be used where p ∈ Z + is the number of monomial terms of u and ξ in the model, c i ∈ R are the co-efficients, b i,j and d i,q are the exponents of the u k−j and ξ k−q factors in the i th monomial term respectively and g r is the linear co-efficient of the r th output term y k−r . Lemma 1: For a model of form (23), and under the assumption that ξ k ∼ N (0, 1) is independent of the input u l for all l ∈ [1, N ], the simulation model, computed in the sense of (4), is given by where P e := {i ∈ [1, p] | d i,q is even ∀q ∈ [1, n ξ ]} and (n − 1)!! := n!
Note that (25) is of the form (17). Define partitions P e := {i ∈ [1, p] | d i,q is even ∀q ∈ [1, n ξ ]} and P o := {i ∈ [1, p] | i / ∈ P e }. This yields Taking the expectation with respect to the noise process ξ, and recognizing that the second term on the right hand side drops out (due to (12)), we get the following which is the desired result. Notice that in the derivation of (24), no approximations were made, and hence the result is exact, as per the definition of the simulation model in (4). However, this will not be the case for the general polynomial NARMAX model, as will be shown in the next Section.

B. General polynomial NARMAX models
We now consider the full polynomial NARMAX model as shown in (1). The procedure to compute a simulation response of a stochastic model proposed in Section IV-A allows us to compute the expectation of the noise terms in (1) systematically. However, the proposed method does not deal with the non-linear dependence of the output signal y on the noise process ξ. Although {u k } is a known sequence and ξ k is distributed normally, the random variable y k is typically not distributed normally. Hence, model terms involving random variables y k cannot be equivalently represented in terms of Hermite polynomials of y k (the inner-product in (9) no longer corresponds to the expectation with respect to the distribution of y k ). Consequently, in order to compute the conditional expectation of (1), one must recursively eliminate all y k terms, which yields a non-linear infinite impulse response (NIIR) in terms of the input. Computing an NIIR is intractable because next to infinite time lags, it also contains infinite order polynomial exponents. As a result, just like in the LTI case, an approximation must be made in order to keep the simulation model tractable. Several approximation concepts are proposed in the following example.
Example 2: Consider a system governed by the following equation where, u k ∼ N (µ u , σ u ) and ξ k ∼ N (0, 1). By taking the expectation and recursively replacing the y k−1 term in the model, we get the following set of equations From these equations, it can be seen that realization of the simulation model results in an NIIR. We explore three approximation concepts, the first two are commonly used in the literature, while the third is our proposed approximation.
(i) Ignore the noise terms in (1). This yields (ii) Truncate the NIIR after a certain number of recursive substitutions of the y k terms. We get the simulation responses for one and two recursive substitutions, respectively, followed by truncation of the NIIR. (iii) After a certain number of recursive substitutions, approximate the tail of the truncated NIIR by the past sample of the simulation response. Applying this concept after one recursive replacement yields In Section V it will be shown that the simulation model (32) approximates E[y k ] better than the other models, in terms of the RMS error.
The purpose of approximation concept (iii) is to achieve a compromise between the approximation concept (i) and the NIIR realiztion of the simulation response. Introduce a parameter l ∈ Z to denote the number of successive recursions of the past output terms before the remaining NIIR is approximated with the past simulation response. Denote the resulting simulation approximation as y s,k,l . Observe that the simulation model in (32) would be labelled as y s,k,1 as per the new notation. The parameter l achieves a trade-off between the NIIR realization and the approximation in (29). As l approaches ∞, the approximation y s,k,l approaches the conditional expectation E[y k ]. Furthermore, for l = 0, the approximation y s,k,0 reduces to (29).
We now formalize the proposed approximation method for the polynomial NARMAX model class.  (1) will be re-written as The l recursive substitutions can be represented in the following steps It should be noted that the functions f, f 1 , . . . f l are not identical. However, they are equivalent, in the sense that for given sequences {u k ; k ∈ [0, N ]} and {ξ k ; k ∈ [0, N ]}, these functions produce the same output sequence {y k ; k ∈ [n, N ]}, wheren = max(n u + l, n y + l, n ξ + l). Moreover, since f is a polynomial function, so are f 1 , . . . , f l . The parameterized approximation can now be defined as follows.

Definition 2:
For the polynomial NARMAX model in (33), the l−approximate simulation model is defined as To compute the l−approximate simulation model, the following alternate representaiton of (34) is used where a i,r is the exponent of y k−r in the i th term.
Theorem 1: The l−approximate simulation model for the polynomial NARMAX model in (1), under the assumption that ξ is independent of u and that ξ k ∼ N (0, 1), is given by Proof: The proof makes use of the representation in (36), the definition of l−approximate simulation, and the use of Hermite polynomials of the noise process as in Section IV-A. The proof follows the same line of reasoning as in the proof of Lemma 1. and is thus omitted here.
Remark 1: It should be noted that the proof of Lemma 1 critically rests on the assumption that the noise terms in (1) are polynomial, and that ξ k ∼ N (0, 1). While Lemma 1 and Theorem 1 are derived for models of the form (23) and (36) respectively, the input terms u and output terms y were not explicitly required to be in a polynomial form. Hence, the derived results can easily be extrapolated to compute simulation models from prediction models with arbitrary non-linear functions on the input and output terms, as long as the noise terms appear in a polynomial form.
Remark 2: While we assumed that ξ k ∼ N (0, 1), the results can be extended to a non-standard Gaussian distribution (see (5)). Furthermore, the results can also be exteded to any exponential family distribution by suitably changing the probability measure and the Hilbert space of random functions in Section III.

V. NUMERICAL ILLUSTRATION
A simulation study was carried out to verify the simulation models derived in Examples 1 and 2. The system in (19) was simulated with a periodic input u ∼ 3 + N (0, 1) with p = 1024 periods of N = 3000 samples, and with noise disturbance ξ ∼ N (0, 1). Data from the first 5 experiments was discarded to remove any errors due to transients. Assuming ergodicity, the ideal simulation response can be computed empirically as the ensemble averagē   where y (i) is the output corresponding to the i th period of the input. The theoretical simulation of the model is computed using (21). The empirical and the theoretical simulation response is plotted in Fig. 1a. Several realizations of the noisy output are plotted in yellow. The unscaled histogram ofȳ s − y s is depicted in Fig. 1b. The differencē y s − y s appears to be centered at 0. This clearly implies that the simulation model computed by neglecting the noise contributions (see (22)) would be biased by a scalar factor of 1. The difference between the two simulation models is also quantified in terms of the RMS error in Table I. The first column of numbers indicate the RMS error between the empirical mean and the proposed simulation models. The second column indicates ensemble average of the RMS errors between the noisy response of the system and the simulation model approximations. It can be observed that the proposed simulation model performs significantly better in approximating the noisy output than the conventional simulation model that sets the noise contributions to 0. A similar simulation study is carried out for the model in Example 2 (see Equation (28)). The model was simulated with periodic input u k ∼ N (0, 1) with p periods of N samples each, and noise ξ k ∼ N (0, 1). Again, the simulation response was computed empirically using p − 5 periods of excitation. The histograms of the differences between the four simulation models (29-32) and the empirical simulation y s are depicted in Fig. 2. Additionally, the RMS errors are presented in Table II. It can be verified, both visually and numerically, that y (4) s approximates the simulation of the system better than the other approximations.

VI. CONCLUSIONS
When computing a simulation model from a prediction model that belongs to the NARMAX class, care must be taken to avoid bias. To that extent, we proposed an alter- native representation that allows to compute a simulation model from a polynomial NARMAX model conveniently.
In order to curtail infinite order simulation models, we also proposed a parameterized method to approximate the simulation model. The proposed approximation method yields the true, infinite-order simulation model as the approximation parameter goes to ∞. Simulation examples were used to demonstrate that the approximation method yields a simulation model with a smaller bias compared to the other common approaches used for computing simulation models.