Novel Insights into Estimation of Bilinear Time Series Models with Exponential and Symmetric Coefficients

Abstract

This paper focuses on the estimation and simulation of a specific subset of bilinear time series models characterized by dynamic exponential coefficients. Employing an exponential framework, we delve into the implications of the exponential function for our estimation process. Our primary aim is to estimate the coefficients of the proposed model using exponential coefficients derived from time-varying parameters. Through this investigation, our goal is to shed light on the asymptotic behaviors of the estimators and scrutinize their existence and probabilistic traits, drawing upon the foundational theorem established by Klimko and Nilsen. The least squares approach is pivotal in both estimating coefficients and analyzing estimator behavior. Moreover, we present a practical application to underscore the real-world implications of our research. By offering concrete examples of applications and simulations, we endeavor to provide readers with a comprehensive understanding of the implications of our work within the realm of time series analysis, specifically focusing on bilinear models and time-varying exponential coefficients. This multifaceted approach underscores the potential impact and practical relevance of our findings, contributing to the advancement of the field of time series analysis. To discern the symmetry characteristics of the model, we estimate it using coefficients that sum to zero and conduct a brief comparative analysis of two bilinear models.

1. Introduction

Symmetry is a key element in analyzing nonlinear time series models, significantly aiding in the extraction of a broad array of properties. It provides an in-depth conceptual understanding and practical benefits for both modeling and forecasting. The presence of symmetry, whether in model coefficients or structural dynamics, simplifies the models by reducing the number of parameters to estimate, thereby decreasing the complexity and enhancing the interpretability. Furthermore, symmetry finds application in error correction and the detection of anomalies. Deviations from expected symmetric patterns in time series can highlight errors, outliers, or extraordinary occurrences, necessitating additional analysis or corrective measures. In light of the constraints posed by classical linear time series models in comprehensively capturing complex phenomena, there has been a notable shift towards the exploration of nonlinear alternatives. Among these, bilinear models of time series have emerged as a focal point, drawing significant interest owing to their remarkable versatility across a spectrum of disciplines. These models have proven their value in diverse fields, including biology, where they elucidate intricate ecological dynamics; medicine, where they aid in understanding nonlinear physiological responses; queueing theory, where they enhance the modeling of service systems; chemistry, where they unravel complex chemical reactions; software reliability, where they predict inter-failure patterns; and signal processing, where they extract valuable information from noisy data streams. This multifaceted applicability underscores the indispensability of nonlinear models, particularly bilinear ones, in tackling the intricacies of real-world phenomena.

This paper builds upon a rich history of research, tracing the evolution of bilinear models. Granger and Andersen’s seminal work in 1978 (as documented in references [1,2,3,4]) marked a significant milestone in the development of these models. Their groundbreaking contributions laid the foundation for an understanding of the intricate dynamics of time series data, introducing a framework capable of capturing the nonlinearity in various phenomena.

Granger and Andersen’s bilinear models offer a unique approach by allowing for interactions between different variables that may evolve. This flexibility proves crucial in addressing the real-world complexities that classical linear models struggle to capture. The framework presented by Granger and Andersen not only provided a theoretical foundation but also opened avenues for practical applications across disciplines.

Subsequent advancements by researchers such as Subba Rao in 1981, Subba Rao and Gabr Tong in 1981, Quinn in 1982, Hannan in 1982, Bhaskara Rao et al. in 1983, Liu and Brockwell in 1988, and Liu in 1990 further expanded the understanding of bilinear models. These contributions delved into aspects such as stationarity, invertibility, and the specific properties of bilinear models, enriching the theoretical landscape and enhancing the applicability of these models in diverse domains.

As a response to the evolving nature of many phenomena, this paper addresses the stationarity conditions of time-varying bilinear models, a concern articulated by Bibi and Oyet in 2002 and 2007. The estimation of bilinear models proves particularly challenging, especially when the coefficients vary over time or assume a functional form dependent on multiple variables.

Structured to provide clarity and depth, Section 2 provides preliminary information on bilinear models, emphasizing relevant theories and introducing the Klimko–Nilsen theorem, establishing the necessary stability conditions. In Section 3, we scrutinize a bilinear model with exponential coefficients, assessing the hypotheses of the Klimko–Nilsen theorem and detailing the approach to estimating model coefficients. In Section 4, with applications and numerical illustrations, including MATLAB-generated graphs, our study concludes with general remarks and comments on the broader implications of our research in the realm of time series analysis involving bilinear models with time-varying exponential coefficients [5,6,7].

2. Bilinear Time Series Models with Time-Varying Coefficients

One of the most thorough analyses of bilinear models with time-varying coefficients in the statistical literature is found in Bib’s work. A subset of mathematical models called bilinear time series models is frequently used to analyze and forecast time series data [8,9]. These models are essential in comprehending the intricate linkages that are present in temporal datasets. These models describe the complex dynamics in time series by seamlessly integrating linear and nonlinear components. The phrase “bilinear” describes the addition of two variables multiplied by one another, which typically represent the system’s input and output. We explore Bib’s groundbreaking work in this section, which has been crucial in helping to understand the behavior of bilinear time series models with time-varying coefficients. We hope to improve further our grasp of the complexity related to this class of models by expanding on Bib’s fundamental discoveries. By carefully analyzing Bib’s approach and results, we hope to identify important details that shed light on the time series analysis field as a whole and provide a basis for the debates and analyses that follow in this study.

Definition 1 

([4]). A general form of a bilinear time series models $\left\{x_t,t\in \mathbb{Z}\right\}$ defined on probability space $(\Omega ,A,P$) with order $(p,q,r,s)$ is represented as follows,

$$x\left(t\right)=\sum _{i=1}^p{u}_i\left(t\right)x(t-i)+\sum _{j=1}^r{v}_j\left(t\right)\epsilon (t-j)+\sum _{k=1}^q\sum _{l=1}^s{\gamma }_{kl}\left(t\right)x(t-k)\epsilon (t-l)+\epsilon \left(t\right),$$

(1)

and noted by $bl(p,r,q,s)$, where $\left(u_i\left(t\right),v_j\left(t\right),{\gamma }_{kl}\left(t\right)\right)$ are the time-varying coefficients. $\epsilon \left(t\right)$ is the white noise part, which is not necessarily identically distributed, with mean zero and variance ${\sigma }^2\left(t\right)$. It takes several forms, such as ARCH, GARCH, or COGARCH in continuous cases or processes [10,11].

There are cases whereby $\epsilon \left(t\right)$ is written as

$$\epsilon \left(t\right)=\zeta \left(t\right)\sigma \left(t\right).$$

(2)

For example, the quantity (${\zeta \left(t\right))}_{t\in \mathbb{Z}}$ is a sequence of independent and identically distributed random variables, where

$$\left\{\begin{array}{c}E\left(\zeta \left(t\right)\right)=E\left({\zeta }^3\left(t\right)\right)=0\\ E\left({\zeta }^2\left(t\right)\right)=E\left({\zeta }^4\left(t\right)\right)=0\end{array}\right.$$

(3)

Definition 2 

([12]). ${\widehat{\alpha }}_N$ is defined as an estimator for α if and only if ${\widehat{\alpha }}_N$ is a solution of

$$arg\; min{q}_N\left(\alpha \right).$$

(4)

where N number of observations, $q_N\left(\alpha \right)$ the penalty function expressed by the following expression

$$q_N\left(\alpha \right)=\frac{1}{N}\sum _{i=1}^N{\epsilon }_{\alpha }^2\left(t\right).$$

(5)

Definition 3 

([11]). The method of least squares is based on Taylor’s formula of the second degree: for any α, a parameter that we wish to estimate, it can be written as follows,

$$q_N\left(\alpha \right)=q_N\left({\alpha }_0\right)+(\alpha -{\alpha }_0)\frac{\partial q_N}{\partial {\alpha }^T}\left({\alpha }_0\right)+\frac{1}{2}(\alpha -{\alpha }_0)\frac{{\partial }^2{q}_N}{\partial {\alpha }^T\partial \alpha }\left({\alpha }^{*}\right){(\alpha -{\alpha }_0)}^T.$$

(6)

where T represents the transpose of a matrix; ${\alpha }^{*}$ is an intermediary point between α and ${\alpha }_0$, where $∥(\alpha -{\alpha }_0)∥\le \delta ,$$\delta >0.$

Definition 4. 

We define the orthogonal projection and, denoted by $p_{t/t-1}$, the following difference:

$$\epsilon \left(t\right)=x\left(t\right)-p_{t/t-1}.$$

(7)

Example 1. 

We have the following sample of models with constant coefficients:

$$x\left(t\right)={\gamma }_{0}x(t-p)+{\gamma }_{1}x(t-p)\epsilon (t-1)+\epsilon \left(t\right).$$

(8)

We revisit the renowned theorem of Klimko and Nilsen, which addresses the existence and requisite conditions for estimators. This is the correct view in many models, such as bilinear models, in cases in which the parameters of the models are pure constants. The theorem is articulated as follows.

Theorem 1. 

Let $\left\{x\left(t\right),t\in \mathbb{Z}\right\}$ be a stable process generated by (8). For example, $p_{t/t-1}$ is almost surely twice continuous in an open subset Ω, containing the true value ${\alpha }^0$=$({\gamma }_0,{\gamma }_1)$ of vector α, and let $K,$$K_0$ be constants. The following four assumptions should be satisfied.

$\left(A1\right)$ $E_{{\alpha }^0}{\left\{\frac{{\partial }^2\epsilon \left(t\right)}{\partial {\gamma }_i^2}\right\}}^4\le K,$ $i=0,1.$

$\left(A2\right)$ $E_{{\alpha }^0}{\left\{\frac{{\partial }^2\epsilon \left(t\right)}{\partial {\gamma }_i\partial {\gamma }_j}-E_{{\alpha }^0}\left(\frac{{\partial }^2\epsilon \left(t\right)}{\partial {\gamma }_i\partial {\gamma }_j}\mid F_t\right)\right\}}^2\le K_0,$ $i,j=0,1.$

$\left(A3\right)$$\frac{1}{N}{\sum }_{i=1}^N{E}_{{\alpha }^0}\left\{\left(\frac{{\partial }^2\epsilon \left(t\right)}{\partial {\gamma }_i\partial {\gamma }_j}\mid F_t\right)\right\}$ converge for matrix M, and this matrix is strictly positive for constants.

$\left(A4\right)$$lim\; sup\frac{1}{N\lambda }{\sum }_{i=1}^N∥E_{\alpha ={\alpha }^0}\left(\frac{{\partial }^2\epsilon \left(t\right)}{\partial {\gamma }_i\partial {\gamma }_j}\right)-E_{\alpha ={\alpha }^{*}}\left(\frac{{\partial }^2\epsilon \left(t\right)}{\partial {\gamma }_i\partial {\gamma }_j}\right)∥<\infty ,$ where $∥{\alpha }^0-{\alpha }^{*}∥<\lambda .$

Then, there exists an estimator $\widehat{\alpha },$ where this estimator shows asymptotic behavior such that

$$\widehat{\alpha }\to {\alpha }^{0}whereN\to \infty .$$

Moreover,

$$\left\{\underset{N\to \infty }{lim}\frac{1}{N}\sum E_{\widehat{\alpha }}\left(\frac{{\partial }^2\epsilon \left(t\right)}{\partial \alpha }\frac{{\partial }^2\epsilon \left(t\right)}{\partial {\alpha }^T}\right)-\sum E_{{\alpha }^0}\left(\frac{{\partial }^2\epsilon \left(t\right)}{\partial \alpha }\frac{{\partial }^2\epsilon \left(t\right)}{\partial {\alpha }^T}\right)\right\}\to 0.$$

Proof. 

See ([13]). □

2.1. Bilinear Model Known by the Symmetry of Its Coefficients

This section is devoted to studying and estimating a type of bilinear model known by the symmetry of its coefficients, i.e., with a sum of 0 in its coefficients. In the field of time series, it plays an exceptional role in prediction and estimation and avoids the repetition problem for simulation by extracting a property field according to symmetry. Let the bilinear model have the following symmetric coefficients:

$$x\left(t\right)=ax(t=v)ϵ(t-1)-ax(t-v)+ϵ\left(t\right),v>0.$$

(9)

We will estimate the model coefficients, where the A value is assigned to $a=0.09$ (

Table 1. Estimation of the model coefficients.

	Number of Simulations	$\mathit{s}=250,\mathit{v}=3$
Size (N)	True Values$(\mathit{a},-\mathit{a})$	Estimated Values$(\widehat{\mathit{a}},-\widehat{\mathit{a}})$
50	(0.09, −0.09)	(0.0772, −0.0825)
100	(0.09, −0.09)	(0.0812, −0.0967)
250	(0.09, −0.09)	(0.0887, −0.0792)
500	(0.09, −0.09)	(0.0892, −0.0913)
1000	(0.09, −0.09)	(0.0902, −0.0898)

).

Remark 1. 

In this type of bilinear model, characterized by its symmetric nature, we observe that an increase in the sample size leads to an enhancement in the estimation accuracy as shown in Figure 1. It becomes evident that the estimated values tend to converge closely toward their true counterparts as the sample size grows. Furthermore, as delineated by the proposed model curve, this symmetric behavior becomes apparent.

The symmetry across the x-axis is almost clear, except for some disturbances, which shows that the symmetric coefficients impart this property to the curves.

2.2. Time-Varying Coefficient Construction

For the construction of time-varying coefficients, let the following sample of bilinear models with coefficients be defined with order $bl(p,0,q,1)$

$$x\left(t\right)=u_p(t,a)x(t-p)+{\gamma }_q(t,b)\epsilon (t-1)x(t-q)+\epsilon \left(t\right).$$

(10)

$u_p(t,a)$ and ${\gamma }_q(t,b)$ are two defined functions from ${\mathbb{R}}^n$, ${\mathbb{R}}^m$, respectively, and take their values in $\mathbb{R}.$$x\left(t\right)$ is not related to the previous result, and we assume that the two functions satisfy the general condition for stability

$$u_p^2(t,a)+{\sigma }^2\left(t\right){\gamma }_q^2(t,b)<0.$$

(11)

where $a=(a_1,a_2,\dots ,a_n)$ and $b=(b_1,b_2,\dots ,b_n)$ are two vectors, according to the reference works ([1]), used to determine the recurring expression of models to prove the convergence of series. We use it in our estimation. In situations in which $p=q,$ it is possible to write

$$x\left(t\right)=\sum _{j=1}^{\left[\frac{t}{p}\right]}\prod _{i=1}^{j-1}\left\{u_p(t-ip,a)+{\gamma }_p(t-ip,b)\epsilon (t-ip-1)\right\}\epsilon (t-jp)+\epsilon \left(t\right).$$

(12)

The ▵ subset included in ${\mathbb{R}}^{n+m}$ contains the true value ${\alpha }^0=(a,b)=(a_0,a_1,\dots ,a_n,$$b_0,b_1,\dots ,b_m)$ of $\alpha .$

Then, we can extract the following formula:

$$\begin{array}{ccc}\hfill \epsilon \left(t\right)& =& x\left(t\right)-u_p(t,a)x(t-p)-\sum _{j=1}^{t-1}{(-1)}^j\left\{\prod _{i=1}^{j-1}{\gamma }_p(t-i,b)\right\}\left\{\prod _{i=1}^{j-1}x(t-i-p)\right\}\hfill \\ & & \times \left(x(t-j)-u_p(t,a)x(t-j-p)\right).\hfill \end{array}$$

(13)

$\left[c\right]$ denotes the integer part of c. In the least squares approach, we are obliged to specify the F $\sigma -$field of events. Thus, we consider the observations of the model $\left(x\right(1),x(2),\dots ,x(t\left)\right)$.

3. Estimation Results

The outcome of least squares estimation is the set of parameter values that minimizes the sum of the squared differences between the observed and predicted values in a regression analysis. This estimation method is commonly used in statistics to find the best-fitting line or curve for a given set of data points. If we have specific data or a regression problem, we typically perform least squares estimation using the relevant mathematical formulas or software tools to obtain the specific numerical values for the coefficients. We use the least squares method to prove our theorem in the case of time-varying coefficients, and we use the likelihood method to estimate the coefficients of a sample of bilinear models with exponential coefficients.

3.1. Klimko–Nilsen Theorem in Time-Varying Coefficients

Now, we can prove the Klimko and Nilsen theorem in the context of time-varying constants. The proof of this theorem is based on bounding $p_{t/t-1}$. Thus, in this case, to demonstrate that $\left|p_{t/t-1}\right|$ is bounded, we simply guarantee that $E\left(\left|p_{t/t-1}\right|\right)\le K$, where K is a positive constant, and we have $M=max\left({\epsilon }_{t-pj}\right)$. Then,

$$\begin{array}{ccc}\hfill E\left(\left|p_{t/t-1}\right|\right)& \le & M\sum _{j=1}^{\left[\frac{t}{p}\right]}E\left|\prod _{i=1}^{j-1}E{u}_p(t-ip,a)+{\gamma }_p(t-ip,b)\epsilon (t-ip-1))\right|\hfill \\ & \le & M\sum _{j=1}^{\left[\frac{t}{p}\right]}\prod _{i=1}^{j-1}E\left|u_p(t-ip,a)+{\gamma }_p(t-ip,b)\epsilon (t-ip-1)\right|.\hfill \end{array}$$

Using the Cauchy–Schwarz inequality, we have

$$\begin{array}{ccc}\hfill \Gamma & =& E\left|u_p(t-ip,a)+{\gamma }_p(t-ip,b)\epsilon (t-ip-1)\right|\hfill \\ & \le & {\left\{E{\left|u_p(t-ip,a)+{\gamma }_p(t-ip,b)\epsilon (t-ip-1)\right|}^2\right\}}^{0.5}\hfill \\ & \le & E{\left\{u_p^2(t-ip,a)+{\gamma }_p^2(t-ip,b){\epsilon }^2(t-ip-1)\right\}}^{0.5}.\hfill \end{array}$$

where

$${\left\{u_p^2(t-ip,a)+{\gamma }_p^2(t-ip,b){\epsilon }^2(t-ip-1)\right\}}^{0.5}\le r<1.$$

(14)

Then,

$$\begin{array}{ccc}\hfill E\left(\left|p_{t/t-1}\right|\right)& \le & M\sum _{j=1}^{\infty }\prod _{i=1}^{j-1}r\hfill \\ & \le & \sum _{j=1}^{\infty }{r}^{j-1}.\hfill \end{array}$$

Since the last series is a converging geometric series, we have that $E\left(\left|p_{t/t-1}\right|\right)<K.$ Meanwhile, we have the following derivatives:

$$\frac{\partial {\epsilon }^2\left(t\right)}{\partial \alpha }=-2\epsilon \left(t\right)\frac{\partial p_{t/t-1}}{\partial \alpha }.$$

(15)

It is easy to show that $\frac{\partial p_{t/t-1}}{\partial \alpha }$ is bounded with the repetition of the Cauchy–Schwarz inequality. Then,

$$\begin{array}{ccc}\hfill E\left\{2\epsilon \left(\alpha \right)\frac{\partial p_{t/t-1}}{\partial \alpha }\left(\alpha \right)\right\}& \le & 8K^{2}E\left({\epsilon }^4\left(\alpha \right)\right)\hfill \\ & \le & 8K^{2}M.\hfill \end{array}$$

However, in assumption A2, we have

$$E\left(\frac{\partial {\epsilon }^2\left(\alpha \right)}{\partial {\alpha }_i\partial {\alpha }_j}\right)=2\left\{\frac{\partial p_{t/t-1}}{\partial {\alpha }_i}\left(\alpha \right)\frac{\partial p_{t/t-1}}{\partial {\alpha }_j}\left(\alpha \right)-\epsilon \left(\alpha \right)\frac{{\partial }^2{p}_{t/t-1}}{\partial {\alpha }_i\partial {\alpha }_j}\left(\alpha \right)\right\}.$$

(16)

Then,

$$\begin{array}{ccc}\hfill {\Pi }_1& =& E_{{\alpha }^0}{\left\{\frac{{\partial }^2\epsilon \left(\alpha \right)}{\partial {\alpha }_i\partial {\alpha }_j}-E_{{\alpha }^0}\left(\frac{{\partial }^2\epsilon \left(\alpha \right)}{\partial {\alpha }_i\partial {\alpha }_j}\mid F_t\right)\right\}}^2\hfill \\ & =& E_{{\alpha }^0}{\left\{\frac{\partial p_{t/t-1}}{\partial {\alpha }_i}\left(\alpha \right)\frac{\partial p_{t/t-1}}{\partial {\alpha }_j}\left(\alpha \right)-\epsilon \left(\alpha \right)\frac{{\partial }^2{p}_{t/t-1}}{\partial {\alpha }_i\partial {\alpha }_j}\left(\alpha \right)\right\}}^2\hfill \\ & =& 4E_{{\alpha }^0}{\left\{\frac{\partial p_{t/t-1}}{\partial {\alpha }_i}\left(\alpha \right)\frac{\partial p_{t/t-1}}{\partial {\alpha }_j}\left(\alpha \right)\right\}}^2+16E\left({\epsilon }^2\left(\alpha \right)\right)E{\left\{\frac{{\partial }^2{p}_{t/t-1}}{\partial {\alpha }_i\partial {\alpha }_j}\right\}}^2\hfill \end{array}$$

(17)

which gives ${\Pi }_1<C.$ In A3, we have

$$\begin{array}{ccc}\hfill {\Pi }_2& =& \frac{1}{N}\sum _{i=1}^N{E}_{{\alpha }^0}\left\{\left(\frac{{\partial }^2\epsilon \left(t\right)}{{\alpha }_i{\alpha }_j}\mid F_t\right)\right\}\hfill \\ & =& \frac{1}{N}\sum _{i=1}^N{E}_{{\alpha }^0}\left\{\frac{\partial p_{t/t-1}}{\partial {\alpha }_i}\left(\alpha \right)\frac{\partial p_{t/t-1}}{\partial {\alpha }_j}\left(\alpha \right)-\epsilon \left(\alpha \right)\frac{{\partial }^2{p}_{t/t-1}}{\partial {\alpha }_i\partial {\alpha }_j}\left(\alpha \right)\right\}.\hfill \end{array}$$

We can write

$$R_{ij}=\frac{1}{N}\sum _{i=1}^N{\left\{E_{{\alpha }^0}{\left(\frac{\partial p_{t/t-1}}{\partial {\alpha }_i}\left(\alpha \right)\frac{\partial p_{t/t-1}}{\partial {\alpha }_j}\left(\alpha \right)\right)}^2\right\}}^{0.5}\le C_0.$$

It is clear that the Cauchy–Schwarz inequality makes the expression bounded $R_{ij}\le C_0,$ where $C_0$ is constant.

The proof of A5 comes from

$${\Pi }_2=E_{{\alpha }^0}{\left\{\frac{{\partial }^2\epsilon \left(\alpha \right)}{\partial {\alpha }_i}-\frac{{\partial }^2\epsilon \left(\alpha \right)}{\partial {\alpha }_j}\right\}}^2\ge 0.$$

(18)

We find

$${\Pi }_2=E_{{\alpha }^0}{\left\{\frac{{\partial }^2\epsilon \left(\alpha \right)}{\partial {\alpha }_i}\right\}}^2+E_{{\alpha }^0}{\left\{\frac{{\partial }^2\epsilon \left(\alpha \right)}{\partial {\alpha }_j}\right\}}^2-2E_{{\alpha }^0}\left\{\frac{{\partial }^2\epsilon \left(\alpha \right)}{\partial {\alpha }_i}\frac{{\partial }^2\epsilon \left(\alpha \right)}{\partial {\alpha }_j}\right\}.$$

(19)

Then, the relation emerges in the following form:

$$\begin{array}{ccc}\hfill 2E_{{\alpha }^0}\left\{\frac{{\partial }^2\epsilon \left(\alpha \right)}{\partial {\alpha }_i}\frac{{\partial }^2\epsilon \left(\alpha \right)}{\partial {\alpha }_j}\right\}& \le & E_{{\alpha }^0}{\left\{\frac{{\partial }^2\epsilon \left(\alpha \right)}{\partial {\alpha }_i}\right\}}^2+E_{{\alpha }^0}{\left\{\frac{{\partial }^2\epsilon \left(\alpha \right)}{\partial {\alpha }_j}\right\}}^2\hfill \\ & \le & E_{{\alpha }^0}{\left\{{\left(\frac{{\partial }^2\epsilon \left(\alpha \right)}{\partial {\alpha }_i}\right)}^4\right\}}^{0.5}+E_{{\alpha }^0}{\left\{{\left(\frac{{\partial }^2\epsilon \left(\alpha \right)}{\partial {\alpha }_j}\right)}^4\right\}}^{0.5}.\hfill \end{array}$$

In other terms, we have $E_{{\alpha }^0}$

$${\left(\frac{{\partial }^2\epsilon \left(\alpha \right)}{\partial {\alpha }_i}\right)}^4\le C_2.$$

(20)

Given that $C_2$ is a positive constant, the use of $A1$ leads us to the following conclusion:

$$E_{{\alpha }^0}\left\{\frac{{\partial }^2\epsilon \left(\alpha \right)}{\partial {\alpha }_i}\frac{{\partial }^2\epsilon \left(\alpha \right)}{\partial {\alpha }_j}\right\}\le \frac{1}{2}\sqrt{C_2}\sqrt{C_2}<\infty .$$

(21)

3.2. Likelihood Estimation Approach

The likelihood estimation method offers a consistent method to solve parameter estimation problems. This implies that maximum likelihood estimates are adaptable to a wide range of estimation scenarios. For instance, they can be utilized in reliability analysis for censored data across different censoring models. In this subsection, we estimate the coefficients of a sample of bilinear models with the likelihood method. These coefficients take an exponential form defined by the formula

$$\left\{\begin{array}{c}x\left(t\right)=e^{-\sqrt[p]{{\phi }_a\left(t\right)}}x(t-p)+e^{-\sqrt[p]{{\phi }_b\left(t\right)}}x(t-p)\epsilon (t-1)+\epsilon (t-1)\\ \epsilon \left(t\right)\mathrm{Gassian}\mathrm{variable},p\ge 1\end{array}\right..$$

(22)

This model type satisfies the requisite condition of stationarity, where

$${\left(e^{-\sqrt[p]{{\phi }_a\left(t\right)}}\right)}^2+{\sigma }^2\left(t\right)\left(e^{-\sqrt[p]{{\phi }_b\left(t\right)}}\right)<1.$$

(23)

Time-varying coefficients can be constructed in the form $a=(a_1,a_2)$ and $b=(b_1,b_2).$ Then, the estimated parameter is $\alpha =(a,b)\in {\mathbb{R}}^4$

$$\left\{\begin{array}{c}{\phi }_a\left(t\right)=\frac{{(-1)}^t+1}{2}{a}_1+\frac{-{(-1)}^t+1}{2}{a}_2\\ {\phi }_b\left(t\right)=\frac{{(-1)}^t+1}{2}{b}_1+\frac{-{(-1)}^t+1}{2}{b}_2\end{array}\right..$$

(24)

It is noticed that this function provides alternative coefficients; depending on expression (13), the expression of ${\epsilon }_a\left(t\right)$ will be

$$\begin{array}{ccc}\hfill {\epsilon }_a\left(t\right)& =& x\left(t\right)-e^{-\sqrt[p]{{\phi }_a\left(t\right)}}x(t-p)-\sum _{j=1}^{t-1}{(-1)}^j\left\{\prod _{i=1}^{j-1}{e}^{-\sqrt[p]{{\phi }_b\left(t\right)}}\right\}\left\{\prod _{i=1}^{j-1}x(t-i-p)\right\}\hfill \\ & & \times \left(x(t-j)-e^{-\sqrt[p]{{\phi }_a\left(t\right)}}x(t-j-p)\right).\hfill \end{array}$$

(25)

Subsequently, by substituting the coefficients, the relationship is written as

$$\begin{array}{ccc}\hfill {\epsilon }_{\alpha }\left(t\right)& =& x\left(t\right)-e^{-\sqrt[p]{a_i}}x(t-p)-\sum _{j=1}^{t-1}{(-1)}^j\left\{\prod _{i=1}^{j-1}{e}^{-\sqrt[p]{b_i}}\right\}\left\{\prod _{i=1}^{j-1}x(t-i-p)\right\}\hfill \\ & & \times \left(x(t-j)-e^{-\sqrt[p]{a_i}}x(t-j-p)\right).\hfill \end{array}$$

(26)

If we use $\sigma \left(t\right)=\sigma$ (the variance is fixed with a constant), we define the penalty function

$$f(x\left(t\right),\alpha )=\frac{1}{\sigma \sqrt{2\pi }}exp(\left(-\frac{{\epsilon }_{\alpha }^2\left(t\right)}{2{\sigma }^2}\right)$$

(27)

To estimate the coefficients of the model, we seek solutions $\alpha$ to minimize the product

$$L\left(\alpha \right)=\prod _{i=1}^{N}f(x\left(t\right),\alpha )={\left(2\pi \right)}^{\frac{-N}{2}}{\sigma }^{-N}exp\left(-\frac{1}{2{\sigma }^2}\sum _{i=1}^N{\epsilon }_{\alpha }^2\left(t\right)\right).$$

(28)

The likelihood function of logarithm $ln\left(L\right(\alpha \left)\right)$ produces

$$\begin{array}{ccc}\hfill \Phi \left(\alpha \right)& =& ln\left(L\right(\alpha \left)\right)\hfill \\ & =& \frac{-N}{2}ln\left(2\pi \right)-Nln\left(\sigma \right)-\frac{1}{2{\sigma }^2}\sum _{i=1}^N\frac{{\epsilon }_{\alpha }^2\left(t\right)}{2{\sigma }^2}\hfill \\ & =& \frac{-N}{2}ln\left(2\pi {\sigma }^2\right)-\frac{1}{2{\sigma }^2}\sum _{i=1}^N\frac{{\epsilon }_{\alpha }^2\left(t\right)}{2{\sigma }^2}\hfill \end{array}$$

In this case, we present the technical derivation for model (22)

$$\begin{array}{ccc}\hfill \frac{\partial \Phi \left(\alpha \right)}{\partial \alpha }& =& -\frac{1}{2{\sigma }^2}\sum _{i=1}^N\frac{\partial {\epsilon }_{\alpha }^2\left(t\right)}{\partial \alpha }\hfill \\ & =& -\frac{1}{{\sigma }^2}\sum _{i=1}^N\frac{\partial {\epsilon }_{\alpha }\left(t\right)}{\partial \alpha }{\epsilon }_{\alpha }\left(t\right).\hfill \end{array}$$

For example,

$$\begin{array}{ccc}\hfill \frac{\partial {\epsilon }_{\alpha }\left(t\right)}{\partial a_i}& =& \frac{\partial }{\partial \alpha }\left\{x\left(t\right)-e^{-\sqrt[p]{a_i}}x(t-p)-\sum _{j=1}^{t-1}{(-1)}^j\prod _{i=1}^{j-1}{e}^{-\sqrt[p]{b_i}}\prod _{i=1}^{j-1}x(t-i-p)\right\}\hfill \\ & & \times \left(x(t-j)-e^{-\sqrt[p]{a_i}}x(t-j-p)\right).\hfill \end{array}$$

In this case, we present the derivation technique for the model (5), so, in a situation in which $\alpha =(a,b)\in I\times I\subset {\mathbb{R}}^{*+}\times {\mathbb{R}}^{*+}$

$${\displaystyle \frac{\partial \Phi \left(\alpha \right)}{\partial a}=-\frac{1}{2{\sigma }^2}\sum _{t=1}^N\frac{\partial {\epsilon }_t^2\left(\alpha \right)}{\partial a}=-\frac{1}{{\sigma }^2}\sum _{t=1}^N{\epsilon }_t\left(\alpha \right)\frac{\partial {\epsilon }_t\left(\alpha \right)}{\partial a}.}$$

(29)

For example,

$$\begin{array}{cc}\hfill \frac{\partial {\epsilon }_t\left(\alpha \right)}{\partial a}& =\frac{\partial }{\partial a}\left[{\displaystyle x_t-e^{-\left(\sqrt[p]{a}\right)}{x}_{t-p}-\sum _{j=1}^{t-1}{(-1)}^{j-1}{e}^{-(j-1)\left(\sqrt[p]{b}\right)}\left\{\prod _{i=0}^{j-1}{x}_{t-i-p}\right\}\left\{x_{t-j}-e^{-\left(\sqrt[p]{a}\right)}{x}_{t-j-p}\right\}}\right]\hfill \\ \hfill & {\displaystyle =a^{\frac{1}{p}-1}{e}^{-\left(\sqrt[p]{a}\right)}-\sum _{j=1}^{t-1}{(-1)}^{j-1}{e}^{-(j-1)\left(\sqrt[p]{b}\right)}\left\{\prod _{i=0}^{j-1}{x}_{t-i-p}\right\}\left(a^{\frac{1}{p}-1}{e}^{-\sqrt[p]{a}}{x}_{t-j-p}\right).}\hfill \end{array}$$

which shows that

$${\displaystyle \frac{\partial \Phi \left(\alpha \right)}{\partial b}=-\frac{1}{2{\sigma }^2}\sum _{t=1}^N\frac{\partial {\epsilon }_t^2\left(\alpha \right)}{\partial b}=-\frac{1}{{\sigma }^2}\sum _{t=1}^N{\epsilon }_t\left(\alpha \right)\frac{\partial {\epsilon }_t\left(\alpha \right)}{\partial b}.}$$

(30)

Furthermore,

$$\begin{array}{cc}\hfill \frac{\partial {\epsilon }_t\left(\alpha \right)}{\partial b}& =\frac{\partial }{\partial b}\left[{\displaystyle x_t-e^{-\left(\sqrt[p]{a}\right)}{x}_{t-p}-\sum _{j=1}^{t-1}{(-1)}^{j-1}{e}^{-(j-1)\left(\sqrt[p]{b}\right)}\left\{\prod _{i=0}^{j-1}{x}_{t-i-p}\right\}(x_{t-j}-e^{-\left(\sqrt[p]{a}\right)}{x}_{t-j-p})}\right]\hfill \\ \hfill & {\displaystyle =-\sum _{j=1}^{t-1}{(-1)}^{j-1}\frac{\partial }{\partial b}{e}^{-(j-1)\left(\sqrt[p]{b}\right)}\left\{\prod _{i=0}^{j-1}{x}_{t-i-p}\right\}(x_{t-j}-e^{-\left(\sqrt[p]{a}\right)}{x}_{t-j-p})}\hfill \\ \hfill & {\displaystyle =\sum _{j=1}^{t-1}{(-1)}^{j-1}\left\{(1-j)b^{\frac{1}{p}-1}{e}^{(1-j)\left(\sqrt[p]{b}\right)}\right\}\left\{\prod _{i=0}^{j-1}{x}_{t-i-p}\right\}(e^{-\left(\sqrt[p]{a}\right)}{x}_{t-j-p}-x_{t-j}).}\hfill \end{array}$$

where

$$\frac{\partial }{\partial b}{e}^{-(j-1)\left(\sqrt[p]{b}\right)}=(1-j)b^{\frac{1}{p}-1}{e}^{(1-j)\left(\sqrt[p]{b}\right)}.$$

(31)

The following derivation gives

$$\begin{array}{ccc}\hfill \frac{{\partial }^2\Phi \left(\alpha \right)}{\partial a\partial a}& =& {\displaystyle -\frac{1}{2{\sigma }^2}\sum _{t=1}^N\frac{\partial {\epsilon }_t^2\left(\alpha \right)}{\partial a\partial a}}\hfill \\ & =& {\displaystyle -\frac{1}{{\sigma }^2}\sum _{t=1}^N\left\{{\left(\frac{\partial {\epsilon }_t\left(\alpha \right)}{\partial a}\right)}^2+{\epsilon }_t\left(\alpha \right)\frac{{\partial }^2{\epsilon }_t\left(\alpha \right)}{\partial a\partial a}\right\}.}\hfill \end{array}$$

Expanding on the significance of this method, we can also provide insights into the true value, denoted as ${\Omega }^0$. Subsequently, we seek to solve the following equation to obtain an estimate of this value.

$$\mathbf{L}{(x_i,\Omega )}_{i=1,\dots ,N}+\mathbf{O}{(x_i,\Omega )}_{i=1,\dots ,N}(\widehat{\Omega }-\Omega )=0.$$

(32)

Thus, this equation can be rewritten as

$$\widehat{\Omega }=\Omega -{\mathbf{O}}^{-1}{(x_i,\Omega )}_{i=1,\dots ,N}\mathbf{L}{(x_i,\Omega )}_{i=1,\dots ,N}.$$

(33)

This method is based on approximating the estimated value, where we can apply the Newton–Raphson method; see [10]. The value of this numerical method is the freedom to choose the initial value, but one cannot neglect the condition of the stationarity of the model. We propose $\Omega ={\Omega }^0$; then,

$$\begin{array}{cc}\hfill {\Omega }^1& ={\Omega }^0-{\mathbf{O}}^{-1}{(x_i,\Omega )}_{i=1,\dots ,N}\mathbf{L}{(x_i,\Omega )}_{i=1,\dots ,N}\hfill \\ \hfill {\Omega }^2& ={\Omega }^1-{\mathbf{O}}^{-1}{(x_i,\Omega )}_{i=1,\dots ,N}\mathbf{L}{(x_i,\Omega )}_{i=1,\dots ,N}\hfill \\ \hfill & ⋮\hfill \\ \hfill {\Omega }^q& ={\Omega }^{q-1}-{\mathbf{O}}^{-1}{(x_i,\Omega )}_{i=1,\dots ,N}\mathbf{L}{(x_i,\Omega )}_{i=1,\dots ,N}.\hfill \end{array}$$

The repetition of the iterative values each time can give a better approximation; then, if q tends to infinity, ${\Omega }^q$ will converge to the estimated value $\widehat{\Omega }.$

4. Numerical Simulations

Our model simulation, as depicted in Equation (16), holds the potential to provide insights and address key questions raised in theoretical studies. By implementing our simulation framework, we can illuminate various aspects of the theoretical analysis and offer empirical validation or refinement to existing theoretical propositions. Through rigorous experimentation and analysis, our simulation model can enhance our understanding of complex theoretical concepts. We denote by s the number of simulations, by N the sample size, and by p the order of the proposed model. We apply our simulation to the model defined by expression (22) (

Table 2. Estimation of the coefficients of the bilinear model defined by expression (5).

	Number of simulations $\mathit{s}=200$	$p=1$
Size ($\mathbf{N})$	True values ${\Omega }^0=(a_1,a_2,b_1,b_2)$	Estimated values $\widehat{\Omega }=({\widehat{a}}_1,{\widehat{a}}_2,{\widehat{b}}_1,{\widehat{b}}_2)$
50		(0.0304, 0.3891, 0.1651, 0.1532)
100		(0.0281, 0.4212, 0.1670, 0.1763)
250	$(0.02,0.45,0.12,0.25)$	(0.0217, 0.4348, 0.1608, 0.1880)
500		(0.0144, 0.4403, 0.1634, 0.1914)
1000		(0.0192, 0.4430, 0.1634, 0.2468)
	Number of simulations $\mathbf{s}=$ $\mathbf{500}$	$p=2$
Size ($\mathbf{N})$	True values ${\Omega }^0=(a_1,a_2,b_1,b_2)$	Estimated values $\widehat{\Omega }=({\widehat{a}}_1,{\widehat{a}}_2,{\widehat{b}}_1,{\widehat{b}}_2)$
50		(0.0226, 0.3990, 0.1215, 0.1722)
100		(0.0184, 0.4118, 0.1280, 0.2285)
250	$(0.02,0.45,0.12,0.25)$	(0.0243, 0.4276, 0.1482, 0.1835)
500		(0.0204, 0.4495, 0.1269, 0.2411)
1000		(0.0219, 0.4435, 0.1590, 0.1886)
	Number of simulations $\mathbf{s}=\mathbf{1000}$	$p=3$
Size ($\mathbf{N})$	True values ${\Omega }^0=(a_1,a_2,b_1,b_2)$	Estimated values $\widehat{\Omega }=({\widehat{a}}_1,{\widehat{a}}_2,{\widehat{b}}_1,{\widehat{b}}_2)$
50		(0.0031, 0.3865, 0.1482, 0.1709)
100		(0.0057, 0.4072, 0.1545, 0.1769)
250	$(0.02,0.45,0.12,0.25)$	(0.0092, 0.4287, 0.1524, 0.1811)
500		(0.0176, 0.4391, 0.1570 0.1849)
1000		(0.0187, 0.4436, 0.1569, 0.1884)

).

Comments. We notice that the estimated values tend towards the true values; we also observe that the first case is a better approximation but is overall almost compatible.

Comments and Main Results

Observing the nearly identical nature of the three curves as shown in Figure 2, Figure 3 and Figure 4, it is evident that the set of models with exponential coefficients exhibits asymptotic behavior. Consequently, variations in the parameter p within the initial cases maintain the overarching structure of our model. This characteristic holds significant importance in financial series modeling, especially when dealing with queue changes. The validity of the Klimko–Nilsen theorem’s hypotheses is affirmed. Demonstrating the theorem’s applicability with the proposed models, it is revealed that this particular set of models yields highly accurate coefficient estimates, contingent upon stability conditions. This underscores their efficacy in practical applications requiring robust estimation methodologies [14,15].

5. Conclusions and Perspectives

The estimation of the presented sample of bilinear models with exponential coefficients demonstrates notable efficiency. This efficiency is evident in the preservation of the asymptotic behavior of the estimators. Furthermore, with an increase in the sample size, there is a clear tendency for the estimated values to converge towards the true values.

In the specific simulation scenario with 500 iterations and $p=2$, the best approximation is achieved, resulting in the closest alignment of the estimated values with the true values. Noteworthy compatibility is observed among the three graphs, serving as a clear demonstration of the value inherent in our contribution.

As a continuation of our research, our subsequent paper will delve into an exploration of the behavior of these models, particularly as the parameter p tends towards infinity. This extended investigation aims to provide a comprehensive understanding of the models under increasingly complex conditions, further contributing to the broader knowledge in this field. Regarding the symmetric property of the coefficients in bilinear models, we find that this symmetry also carries over to the curves that they generate. This shows that this property is heritable from the coefficients to the curves.