Learning governing physics from output only measurements

Dynamical processes are generally expressed in terms of higher order differential equations. However, in most of the data-driven machine learning techniques the actual system is not observed through its original space but identified in a projected space. The projected space contains all the states of a system which are directly observable. For instance, the second order dynamical systems are often expressed in terms of its displacement and velocity components. Let the projection be realized by a map $\mathcal{T}:\mathbb{R}^{d}\to\mathbb{R}^{m}$ such that $m<d$. If there exists a mapping $\mathcal{T}$, then any higher order system can be reduced to a set of SDEs of the form:

where ${\bm{f}}({{\bm{X}}_{t}},t)$ is the deterministic dynamics of the underlying process, ${\bm{g}}({{\bm{X}}_{t}},t)$ is the volatility associated with the dynamics arising due to the stochastic external input, and ${\bm{\zeta}}(t)$ is the white noise. A more appropriate representation of the above equation can be derived through first order Itô SDEs which arises naturally in non-linear dynamical systems subjected to stochastic excitation such as earthquake, wind force, wave force etc. [41, 42]. In non-linear stochastic dynamical systems the SDEs provide an splendid approach for expressing the behavior of the underlying dynamical system. The SDEs are generally defined in term of its deterministic drift and the additive stochastic diffusion components where both the components are allowed to depend on time and systems states. Let $\left({\Omega,\mathcal{F},P}\right)$ be the probability space and $\{{\mathcal{F}_{t}},0\leq t\leq T)\}$ be the natural filtration constructed from sub $\sigma$-algebras of $\mathcal{F}$. Under the probability space $\left({\Omega,\mathcal{F},P}\right)$ an $m$-dimensional $n$-factor SDE driven by $n$-dimensional Brownian motion {${\bm{B}}_{j}(t),j=1,\ldots n$} is defined as,

where ${\bm{X}}_{t}\in{\mathbb{R}^{m}}$ denotes the ${{\mathcal{F}_{t}}}$-measurable state vector, ${\bm{f}}\left({{{\bm{X}}_{t}},t}\right):{\mathbb{R}^{m}}\mapsto{\mathbb{R}^{m}}$ is the drift vector, ${\bm{g}}\left({{{\bm{X}}_{t}},t}\right):{\mathbb{R}^{m}}\mapsto{{\mathbb{R}}^{m\times n}}$ is the diffusion matrix and ${{\bm{B}}_{j}}\left(t\right)\in{\mathbb{R}^{n}}$ is the Brownian motion. The white noises are generalized derivative Brownian motion i.e. ${\bm{\zeta}}(t)=\dot{\bm{B}}(t)$; this ensures that there exists a corresponding SDE of Eq. (2). In a compact matrix notation, the SDEs are can be expressed as,

The solution to Eq. (3) can be obtained using various stochastic integration schemes [43, 44, 45, 46]. The discovery of governing physics in terms of an SDE requires the identification of the drift and diffusion components from output measurements. One of the challenges arising in the identification of the SDEs from output data is the non-differentiability of the Brownian motion. In general, the stochastic Brownian motion $\bm{B}(t)$ is not a well defined mathematical function since it is not differentiable everywhere with respect to the process $\bm{X}(t)$. However, its continuity is assumed to exists in mean square sense. Formally, if the interval $s\in[0,t]$ is partitioned into $n$-parts as, ${\mathcal{P}_{n}}(0,t):=\left[{0={s_{0}}<s_{1}<\ldots<{s_{n}}=t}\right]$, then for a random process $w$, the quadratic variation ${Q_{n}}(w,t)=\sum\nolimits_{i}^{n}{{{\left|{w({s_{i}})-w({s_{i-1}})}\right|}^{2}}}$ converges to ${Q_{n}}(w,t)\to Q(w,t)$. From the property of the Brownian motions, it can be found that $Q(w,t)=t$; this states that $B(t)$ has zero finite variation but non-vanishing quadratic variation. On the other hand, the deterministic functions have finite variation but zero quadratic variations. Under these properties, the Kramers-Moyal formula suggests that the drift and diffusion components of an SDE can be extracted from the sample time history in terms of their linear and quadratic variations, respectively. Since the continuity of Brownian motions are mathematically defined in the mean square sense, as a consequence, the diffusion components does not have finite variations but are bounded by their quadratic variations. Thus, the diffusion components, unlike the drifts of an SDE are recoverable only through its covariation terms. The linear and quadratic variations denotes the first and second moments of the sample increments. Leveraging this facts, the drift and diffusion components of an SDE given in Eq. (3) can be expressed in terms of the sample time history as follows (the simplified detailed derivation is presented in A):

where ${{\bm{f}}_{i}}\left({{{\bm{X}}_{t}},t}\right)$ is the $i^{th}$ drift component and ${{\bf{\Gamma}}_{ij}}$ is the $(ij)^{th}$ component of the diffusion covaraince matrix ${\bf{\Gamma}}\in{\mathbb{R}^{n\times n}}:={\bm{g}}(t,{{\bm{X}}_{t}}){\bm{g}}{(t,{{\bm{X}}_{t}})^{T}}$. We assume that the analytical form of the drift and diffusion components can be obtained from measurement sample paths in terms of some basis functions. Let ${\upsilon_{k}}({{\bm{X}}_{t}})$ represents the various linear and non-linear mathematical functions evaluated on the system states, and $N$ is the length of sample path. Also let ${\bf{L}}^{f}\in\mathbb{R}^{N\times K}:=[\upsilon_{1}^{f}({{\bm{X}}_{t}}),\ldots,\upsilon_{K}^{f}({{\bm{X}}_{t}})]$ and ${\bf{L}}^{g}\in\mathbb{R}^{N\times K}:=[\upsilon_{1}^{g}({{\bm{X}}_{t}}),\ldots,\upsilon_{K}^{g}({{\bm{X}}_{t}})]$ are the library of candidate functions constructed from the set $\{\upsilon_{k}({{\bm{X}}_{t}}),k=1,\ldots,K\}$. Here the candidate functions can contain several functional forms such as polynomial, trigonometric, etc. Then, the $i^{th}$ drift component and the ${ij}^{th}$ term of diffusion covariance matrix can be expressed as linear combination of the library functions as, ${f_{i}}({{\bm{X}}_{t}},t)=\sum\nolimits_{k=1}^{K}{{\upsilon_{k}}({{\bm{X}}_{t}}){\theta_{i,k}}}$; $i=1,\ldots,m$ and ${{\bf{\Gamma}}_{ij}}=\sum\nolimits_{k=1}^{K}{{\upsilon_{k}}({{\bf{X}}_{t}}){{\bf{\theta}}_{ij,k}}}$; $i,j=1,2,\ldots m$, briefly,

where $\theta_{i_{k}}^{f}$ and $\theta_{{ij}_{k}}^{g}$ are the weights associated with the $k^{th}$ basis functions of drift and diffusion covariance components, respectively. In general, the libraries ${\bm{L}}^{f}$ and ${\bm{L}}^{g}$ can be different but one can construct a library by taking into account of all the basis functions for discovery of both drift and diffusion terms. In a general setting, the discovery of the drift terms culminates in the following regression problem,

The problem for discovery of diffusion components can also be formulated in the same way as above, thus omitted here for brevity. As a subset of the above regression problems, one can solve the following equations and identify the SDE components independently,

where ${{\bm{Y}}_{i}}$ and ${{\bm{Y}}_{ij}}$ are the target vectors of the sparse regression problem associated with $i^{th}$-drift component and $(ij)^{th}$-diffusion covariance term, respectively, and, ${{\bm{\varepsilon}}_{i}}$ and ${{\bm{\eta}}_{ij}}$ are the corresponding residual error vectors. For the discovery of the $i^{th}$ drift component in Eq. (7a), the target vector is defined as, ${{\bm{Y}}_{i}}={{\Delta t}^{-1}\left[{\left({{X_{i,1}}-{\xi_{i,1}}}\right)\ldots\left({{X_{i,N}}-{\xi_{i,N}}}\right)}\right]^{T}}$, and in the formulation of the sparse regression problem for diffusion components in Eq. (7b), the target vector is set as, ${{\bm{Y}}_{ij}}={{\Delta t}^{-1}\left[{\left({{X_{i,1}}-{\xi_{i,1}}}\right)\left({{X_{j,1}}-{\xi_{j,1}}}\right)\ldots\left({{X_{i,N}}-{\xi_{i,N}}}\right)\left({{X_{j,N}}-{\xi_{j,N}}}\right)}\right]^{T}}$. The weights vectors for $i^{th}$ drift component and ${ij}^{th}$ term of diffusion covariance matrix are given as, ${{\bm{\theta}}_{i}^{f}}={\left[{{\theta_{1_{1}}}},{{\theta_{1_{2}}}},\ldots,{{\theta_{1_{K}}}}\right]^{T}}$, and ${\bm{\theta}}_{ij}^{g}={\left[{\theta_{{ij}_{1}}},{\theta_{{ij}_{2}}},\ldots,{\theta_{{ij}_{K}}}\right]^{T}}$, respectively. The straightforward application of the Gibbs sampling mentioned in the Eqs. (17) $\to$ (22) can be performed to solve the regression problems in Eq. 7. Noting that a $m$-dimensional process has $m$ numbers of drift components it requires $m$-independent regressions to identify all the drifts. Further it can be noted that in a $m$-dimensional process with $m\times n$-diffusion matrix the diffusion covariance matrix ${\bf\Gamma}$ has a dimension $m\times m$. However, due to the symmetry of the covariance matrix, it requires ${{m(m+1)}}/{2}$ numbers of independent regressions to be performed to completely identify the diffusion space.

Since Eq. (7a) and (7b) reflects the same sparse regression problem to be solved, for the brevity they can be represented using the following one-dimensional regression equation:

where ${\bm{Y}}\in\mathbb{R}^{N}$ denotes the $N$-dimensional target vector, ${\bf{L}}$ denotes the library of candidate functions, ${\bm{\theta}}$ is the weight vector, and $\epsilon\in\mathbb{R}^{N}$ is the residual error vector representing the model mismatch error. In the Bayesian inference, given the output measurements the aim is to find the posterior distribution of the weight vector ${\bm{\theta}}$. For estimating the posterior distribution of the weight vector ${\bm{\theta}}$, one can apply the Bayes formula in Eq. (8), which yields,

Modeling the mismatch error $\epsilon$ as i.i.d Gaussian random variable with zero mean and variance $\sigma^{2}$, the likelihood function is written as,

where ${\bf{I}}_{N\times N}$ denotes the ${N\times N}$ identity matrix. The discovery of governing physics demands that the resulting model should be interpretable i.e. the solution of weight vector ${\bm{\theta}}$ should contain only few terms in the final model. In this work, the discontinuous spike and slab (DSS) distributions are used for promotion of sparsity in the solution. The equation discovery using DSS-prior requires the classification of weights of the basis functions into either of the spike and slab component. This is done by introducing a latent indicator variable ${\bm{Z}}=\left[Z_{1},\ldots,Z_{K}\right]$ for each of the component $\theta_{k};k=1,\ldots,K$ of the weight vector ${\bm{\theta}}$. The latent indicator variable $Z_{k}$ is assigned a value 1 if the weight corresponds to the slab component, otherwise, it takes a value 0 when the weight belongs to spike component. It is to be noted that the components of the weight vector that belongs to spike does not contribute to discovery of governing physics. Therefore, let us consider a vector ${\bm{\theta}}_{r}\in\mathbb{R}^{r}:\{r\ll K\}$, composed from the elements of the weight vector ${\bm{\theta}}$ for which $Z_{k}=1$. Denoting ${\bm{\theta}}_{r}$ as the weight vector containing only those variables from ${\bm{\theta}}$ for which $Z_{k}=1$, the DSS-prior is defined as [36, 38],

where the spike and slab distributions are defined as, ${p_{spike}}({\theta_{k}})={\delta_{0}}$ and ${p_{slab}}({{\bm{\theta}}_{r}})=\mathcal{N}\left({{\bf{0}},{\sigma^{2}}{\vartheta_{s}}{{\bf{R}}_{0,r}}}\right)$ with ${{\bf{R}}_{0,r}}={{\bf{I}}_{r\times r}}$. It can be observed that the spike and slab distributions are assumed to be independent. To increase the effectiveness and faster convergence to the sparse solution the noise variance ${\sigma^{2}}$ and the slab variance ${\vartheta_{s}}$ are treated as random variables. The noise variance $\sigma^{2}$ is assigned the Inverse-gamma distribution with the hyperparameters ${\alpha_{\sigma}}$ and ${\beta_{\sigma}}$, the the slab variance $\vartheta_{s}$ is assigned the Inverse-gamma prior with the hyperparameters ${\alpha_{\vartheta}}$ and ${\beta_{\vartheta}}$, the latent variables $Z_{k}$ takes a value from the set $\{0,1\}$, thus each of the element in the latent vector is assigned the Bernoulli prior with common hyperparameter ${p_{0}}$. The hyperparameter ${p_{0}}$ is simulated from the Beta prior with the hyperparameters ${\alpha_{p}}$ and ${\beta_{p}}$. This can be summarized as,

With all the random variables, we construct a bigger Bayesian hierarchical model as shown in Fig. 2, where the hyperparameters ${\alpha_{\vartheta}}$, ${\beta_{\vartheta}}$, ${\alpha_{p}}$, ${\beta_{p}}$, ${\alpha_{\sigma}}$, and ${\beta_{\sigma}}$ are provided as a deterministic constants in the hierarchical model. From the DAG structure in Fig. 2, the joint distribution of the random variables ${\bm{\theta}}$, ${\bm{Z}}$, ${\vartheta_{s}}$, ${\sigma^{2}}$, ${p_{0}}|{\bm{Y}}$ are obtained from the Fig. 2 as,

where $p\left({{\bm{\theta}},{\bm{Z}},{\vartheta_{s}},{\sigma^{2}},{p_{0}}|{\bm{Y}}}\right)$ denotes the joint distribution of the random variables, $p\left({{\bm{Y}}|{\bm{\theta}},{\sigma^{2}}}\right)$ denotes the likelihood function, $p\left({{\bm{\theta}}|{\bm{Z}},{\vartheta_{s}},{\sigma^{2}}}\right)$ is the prior distribution for the weight vector ${\bm{\theta}}$, $p\left({{\bm{Z}}|{p_{0}}}\right)$ is the prior distribution for the latent vector ${\bm{Z}}$, $p\left({{\vartheta_{s}}}\right)$ is the prior distribution for the slab variance ${\vartheta}_{s}$, $p\left({{\sigma^{2}}}\right)$ is the prior distribution for the noise variance, $p\left({{p_{0}}}\right)$ is the prior distribution for the success probability $p_{0}$ and ${p\left({\bm{Y}}\right)}$ is the marginal likelihood. Direct sampling from the joint distribution function is intractable in this case due to the spike and slab distribution. Thus, Gibbs sampling technique is used to draw the random samples from the joint distribution [47]. For the Gibbs sampling, the conditional distributions of the random variables are derived in B. The weights corresponding to the spike distribution does not contribute to the selection of the correct basis functions. Thus, only the weights corresponding to $Z_{k}\neq 0$ i.e. the ${\bm{\theta}}_{r}$ vector is sampled using the Gibbs sampling. Referring to the Eqs. (12), (13), (14), and (15), the sequence of the random variables ${{\bm{\theta}}^{(0)}},{\sigma^{2(0)}},\vartheta_{s}^{(0)},p_{0}^{(0)},{{\bm{Z}}^{(0)}},\ldots,{{\bm{\theta}}^{(1)}},{\sigma^{2(1)}},\vartheta_{s}^{(1)},p_{0}^{(1)},{{\bm{Z}}^{(1)}},\ldots$, using the Gibbs sampling technique is obtained by following steps,

1. 1.

   The weight vector ${\bm{\theta}}_{r}^{(i)}$ is sampled from the Gaussian distribution with mean ${{\bm{\mu}}_{\theta}}$ and variance ${{\bf{\Sigma}}_{\theta}}$ as,

   $${{\bm{\theta}}_{r}^{(i)}}|{\bm{Y}},{\vartheta_{s}^{(i)}},{\sigma^{2(i)}}\sim\mathcal{N}\left({{{\bm{\mu}}_{\theta}^{(i)}},{{\bf{\Sigma}}_{\theta}^{(i)}}}\right),$$

   (17)

   where the mean and covariance is defined as, ${{\bm{\mu}}_{\theta}^{(i)}}={\bf{\Sigma}}_{\theta}^{(i)}{\bf{L}}_{r}^{(i)T}{\bm{Y}}$ and ${{\bf{\Sigma}}_{\theta}^{(i)}}={\sigma^{2(i)}}{\left({{\bf{L}}_{r}^{(i)T}{{\bf{L}}_{r}^{(i)}}+\vartheta_{s}^{{(i)}-1}{\bf{R}}_{0,r}^{{(i)}-1}}\right)^{-1}}$, respectively.

2. 2.

   The latent variable $Z_{k}^{(i+1)}$ is assigned the values from the set, $\left\{0,1\right\}$ by using the Bernoulli distribution as,

   $${Z_{k}^{(i+1)}}|{\bm{Y}},{\vartheta_{s}^{(i)}},{p_{0}^{(i)}}\sim Bern\left({{u_{k}}}\right),$$

   (18)

   where ${u_{k}}=\dfrac{{{p_{0}}}}{{{p_{0}}+\lambda\left({1-{p_{0}}}\right)}}$ and $\lambda=\dfrac{{p\left({{\bm{Y}}|{Z_{k}^{(i)}}=0,{{\bm{Z}}_{-k}^{(i)}},{\vartheta_{s}^{(i)}}}\right)}}{{p\left({{\bm{Y}}|{Z_{k}^{(i)}}=1,{{\bm{Z}}_{-k}^{(i)}},{\vartheta_{s}^{(i)}}}\right)}}$. Here, ${{\bm{Z}}_{-k}^{(i)}}\in\mathbb{R}^{K-1}$ denotes the latent variable vector ${\bm{Z}}$ consisting of all the elements except the $k^{th}$ component. The probability that the $k^{th}$ latent variable $Z_{k}^{(i)}$ takes a value 0 or 1 is estimated as follows,

   $$\begin{array}[]{ll}p\left({{\bm{Y}}|{\bm{Z}}^{(i)},{\vartheta_{s}^{(i)}}}\right)&=\dfrac{{\Gamma\left({{\alpha_{\sigma}}+\dfrac{N}{2}}\right)\beta_{\sigma}^{{\alpha_{\sigma}}}}}{{\Gamma\left({{\alpha_{\sigma}}}\right){{\left({2\pi}\right)}^{\dfrac{N}{2}}}{{\left({{\beta_{\sigma}}+\dfrac{1}{2}{{\bm{Y}}^{T}}{\bm{Y}}}\right)}^{\left({{\alpha_{\sigma}}+\dfrac{N}{2}}\right)}}}};\quad\text{when all $\left\{{Z_{k}^{(i)}}{:k=1,\ldots,K}\right\}=0$}\\ &=\dfrac{{\Gamma\left({{\alpha_{\sigma}}+\dfrac{N}{2}}\right)\beta_{\sigma}^{{\alpha_{\sigma}}}{{\left({\left|{{\bf{R}}_{0,r}^{{(i)}-1}}\right|\left|{{{\bf{\Sigma}}_{\theta}^{(i)}}}\right|}\right)}^{\dfrac{1}{2}}}}}{{\Gamma\left({{\alpha_{\sigma}}}\right){{\left({2\pi}\right)}^{\dfrac{N}{2}}}\vartheta_{s}^{\dfrac{{{h_{z}}}}{2}}{{\left({{\beta_{\sigma}}+\dfrac{1}{2}{{\bm{Y}}^{T}}\left({{{\bf{I}}_{N\times N}}-{{\bf{L}}_{r}^{(i)}}{{\bf{\Sigma}}_{\theta}^{(i)}}{\bf{L}}_{r}^{(i)T}}\right){\bm{Y}}}\right)}^{\left({{\alpha_{\sigma}}+\dfrac{N}{2}}\right)}}}};\quad\text{otherwise.}\end{array}$$

   (19)

3. 3.

   The noise variance ${\sigma^{2(i+1)}}$ is simulated from the Inverse-gamma distribution as,

   $${\sigma^{2(i+1)}}|{\bm{Y}},{\bm{Z}^{(i+1)}},{\vartheta_{s}^{(i)}}\sim IG\left({{\alpha_{\sigma}}+\dfrac{N}{2},{\beta_{\sigma}}+\dfrac{1}{2}\left({{{\bm{Y}}^{T}}{\bm{Y}}-{\bm{\mu}}_{\theta}^{(i)T}{\bf{\Sigma}}_{\theta}^{{(i)}-1}{{\bm{\mu}}_{\theta}^{(i)}}}\right)}\right).$$

   (20)

4. 4.

   The slab variance ${\vartheta_{s}^{(i+1)}}$ is sampled from the Inverse-gamma distribution as,

   $${\vartheta_{s}^{(i+1)}}|{\bm{\theta}^{(i)}},{\bm{Z}}^{(i+1)},{\sigma^{2(i+1)}}\sim IG\left({{\alpha_{\vartheta}}+\dfrac{{{h_{z}}}}{2},{\beta_{\vartheta}}+\dfrac{1}{{2{\sigma^{2}}}}{\bm{\theta}}_{r}^{(i)T}{\bf{R}}_{0,r}^{{(i)}-1}{{\bm{\theta}}_{r}^{(i)}}}\right).$$

   (21)

5. 5.

   The success rate ${p_{0}^{(i+1)}}$ is sampled from the Beta distribution as,

   $${p_{0}^{(i+1)}}|{\bm{Z}}^{(i+1)}\sim Beta\left({{\alpha_{p}}+{h_{z}},{\beta_{p}}+K-{h_{z}}}\right),$$

   (22)

   where ${h_{z}}=\sum\nolimits_{k=1}^{K}{{Z_{k}}^{(i+1)}}$.

6. 6.

   The weight vector ${\bm{\theta}}_{r}^{(i+1)}$ is updated using the step 1.

This MCMC is performed for a total of ${N_{M}}$ simulations out of which initial 1000 samples are discarded as the burn-in samples. Let ${N_{s}}$ denote the number of MCMC required to achieve the stationary distribution after the burn-in samples are discarded. Then the marginal posterior inclusion probability (PIP):= $p\left({{Z_{k}}=1|{\bm{Y}}}\right)$ for each of the $K$ basis functions can be estimated by taking mean over the Gibbs samples for each of the $k^{th}$ latent vector ${Z_{k}}$ [36] as,

The basis functions whose corresponding PIP values are more than 0.5 i.e. $p\left({{Z_{k}}=1|{\bm{Y}}}\right)>0.5$ are included in the final model of discovered equation. The PIP value greater than 0.5 indicates that the selected basis functions are observed more than half of the times in the MCMC simulations. Higher PIP value suggests that in case of unseen scenarios, the corresponding basis functions are highly likely to occur in the data representation target vector. The mean and covariance of the weight vector gives the expected value and standard deviation of the system parameters. The standard deviation provides the confidence interval of the parameters which can be used to design the lower and upper bounds of a random variable in an unseen environment. A pseudo code for the proposed framework, is provided in the Algorithm 1.

The Gibbs sampler is initialised with the following initial values of the hyperparameters: $p_{0}^{(0)}$=0.1, $\vartheta^{(0)}$=10, and $\sigma^{2(0)}$ is set equal to the residual variance from ordinary least-squares regression. The initial vector of binary latent variables ${\bm{Z}}^{(0)}$ is computed by setting ${\bm{Z}}^{(0)}=[Z_{1},\ldots,Z_{K}]$ to zero and then activating the components of ${\bm{Z}}$ that reduce the mean-squared error between the training data and the obtained model from ordinary least-squares. For this purpose a forward followed by a backward search algorithm is devised, where the backward search iterates through the activated components of initial latent vector in similar fashion to forward search. Given all the other parameters, the initial value of $\theta^{(0)}$ is obtained from the Eq. (16). For the commencement of the algorithm, the deterministic prior parameters are set to the following values: $a_{p}$=0.1 and $b_{p}$=1 for the Beta prior on $p_{0}$, $a_{v}$=0.5 and $b_{v}$=0.5 for inverse-Gamma prior on slab variance, and, $a_{\sigma}$=$10^{4}$ and $b_{\sigma}$=$10^{4}$ for inverse-Gamma prior on measurement noise. A Markov chains with 3000 samples is used for Gibbs sampling. The first 1000 samples are discarded as burn-in, and the remaining 2000 samples are used for posterior computation.

A Black–Scholes SDE (formulated as geometric Brownian motion) is a continuous-time stochastic process where the logarithm of the randomly varying quantity follows a Brownian motion. The Black-Scholes SDE is frequently used in stock market for modelling of the evolution of stock price of an underlying asset [48]. The Black–Scholes SDE has the drift $f(t,X_{t})=\lambda X$ and diffusion $g(t,X)=\mu X$ with $\lambda>0$ and $\mu>0$ being the real constants. Towards this, the Black–Scholes SDE is defined as a geometric Brownian motion as follows:

where $B=\{B(t);t\geq 0\}$ is the Brownian motion. Data is generated by simulating Eq. (25) and then corrupting it with Gaussian white noise. Identifying the diffusion term as, $g({\bf{X}}_{t},t)=\mu X(n)$, the variance of the diffusion is obtained: $\Gamma=\mu^{2}X^{2}(n)$. The results of this case study are plotted in the Figs. 3(a) and 4(a). Fig. 3(a) depicts that there are two identified basis function $X$ and $|X|$. However, from the Eq. (25), it is obvious that only one function $X$ should have been identified as the driving function of the Black-Scholes SDE. To understand this discrepancy one needs to consider that the evolution of the random variable $X(t)$ in the Black-Scholes SDE follows Geometric Brownian motion (GBM) and GBMs always assume positive values (e.g. real stock price). Since the characteristics of both the functions $X$ and $|X|$ in this case are same, the algorithm identifies both the functions with almost equal probability. A similar phenomenon is observed in case of the identification of the diffusion term too. Figs. 4(b) depicts almost equal contribution from the functions $X^{2}$ and $X|X|$, however, only the term $X^{2}$ should have been actually identified. Upon considering the corresponding parameter values of the basis functions $X$ and $|X|$ for the drift component, and, $X^{2}$ and $X|X|$ for the diffusion component, the Black-Scholes SDE in Eq. (25) can be easily simulated for any unseen environmental conditions.

To verify that this is not a limitation of the proposed scheme, this numerical demonstration is repeated in the Figs. 3(b) and 4(b) without considering the functions $|X|$ and $X|X|$ in the library. The parameters of the basis functions for the drift and diffusion components are portrayed in Figs. 5(c) and 5(d), respectively. The results clearly show that the proposed scheme is able to identify the basis functions along with their associated parameters $\lambda$ and $\mu$ without any significant error. As a consequence of the above case study, one can ask a question about what basis functions are to be considered in the library. One of the possibilities could be to visualize the time series data and then take scientific and engineering judgements. For example, if the time history data shows no zero crossing then one can opt against considering the functions which shares similar properties (for e.g. $|X|$ and $X$, and $X|X|$ and $X^{2}$). Besides one can also choose to consider due to the fact that the similar terms will have equal probability of occurrence and as a combination will demonstrate the observed phenomenon. However, for future prediction the model will be applicable to process depicting no zero crossing only.

The second order non-linear hardening Duffing-Van der pol oscillator is considered, in this section. The Duffing-Van der pol oscillator draws its physical relevance from the models in flow-induced structural vibration problems. The Duffing-Van der pol oscillator has a cubic dissipating force and it is driven by multiplicative white noise. Towards this, the equation of motion of the undertaken system is expressed as:

where $m$, $c$, and $k$ are the mass, damping, and stiffness parameters of the oscillator. Here, $\alpha$ is a real-valued parameter associated with the cubic non-linearity, $\sigma\geq 0$ is the strength of the multiplicative white noise, and $B=\{B(t);t\geq 0\}$ is the Brownian motion. The time derivative of the Brownian $\dot{B}(t)$ here represents the white noise. With a statespace transformation $X={X_{1}}$, and $\dot{X}={X_{2}}$, where $X_{1}$ and $X_{2}$ represents the displacement and velocity, the corresponding first order Itô-stochastic SDEs for the dynamical system are derived as,