A Three-Stage Nonparametric Kernel-Based Time Series Model Based on Fuzzy Data

Abstract

In this paper, a nonlinear time series model is developed for the case when the underlying time series data are reported by $LR$ fuzzy numbers. To this end, we present a three-stage nonparametric kernel-based estimation procedure for the center as well as the left and right spreads of the unknown nonlinear fuzzy smooth function. In each stage, the nonparametric Nadaraya–Watson estimator is used to evaluate the center and the spreads of the fuzzy smooth function. A hybrid algorithm is proposed to estimate the unknown optimal bandwidths and autoregressive order simultaneously. Various goodness-of-fit measures are utilized for performance assessment of the fuzzy nonlinear kernel-based time series model and for comparative analysis. The practical applicability and superiority of the novel approach in comparison with further fuzzy time series models are demonstrated via a simulation study and some real-life applications.

1. Introduction

The field of time series analysis comprises methods used to analyze the characteristics of a response variable with respect to time. It takes into consideration the fact that observations made over time may have an internal structure (such as autocorrelations, trends, seasonal and/or cyclic variations) that should be accounted for. The main aims of time series analysis are as follows:

• Trend analysis: to identify the underlying pattern or trend in the data over time, such as an upward or downward trend.
• Seasonality analysis: to identify if the data exhibit a repeating pattern over a set period, such as daily, weekly, or yearly.
• Forecasting: to forecast future values using historical data.
• Anomaly detection: to identify any unusual or unexpected observations in the data that deviate from the normal pattern.
• Model selection: to choose an appropriate model to represent the underlying relationships between variables in the data.
• Noise reduction: to remove any unwanted variability or random fluctuations from the data to improve the accuracy of predictions and make the underlying patterns more clear.

These aims can inform decision makers, provide insight into the underlying patterns and relationships in the data, and support the development of data-driven strategies in various fields such as economics, engineering, finance, and more (see, e.g., [1,2,3,4,5,6,7,8]).

Common time series models rely on exact observations and ensure crisp predictions. However, due to various uncertainty factors, it is sometimes preferable to make predictions using imprecise values. For instance, we usually observe imprecise observations in carbon emissions, social benefits and oil reserves, among others [9]. Traditional statistical time series models fail to address prediction problems based on ambiguous or vague information represented by fuzzy data. This shortcoming can be overcome by time series models that use techniques of fuzzy statistics. In general, fuzzy statistics is a branch of statistics that deals with uncertainty and imprecision, e.g., in the data. It includes, for instance, the fields of fuzzy estimation, fuzzy regression, fuzzy clustering, and fuzzy hypothesis testing [10,11,12].

Fuzzy time series models were originally introduced in 1993 [13], and since then they have replaced conventional (crisp) time series approaches when observations are uncertain. When considering fuzzy time series models, the prediction of future values requires three principal steps. In step 1, the exact data are reported. In step 2, through the identification of fuzzy logical relations [14,15], the predictions are transformed into fuzzy quantities. Finally, step 3 provides a defuzzification approach [16,17,18,19,20,21,22] to transform the fuzzy values into crisp ones. The techniques used to identify fuzzy logical relations in step 2 primarily involve fuzzy logical relation groups and matrices [13,23,24,25,26,27,28,29,30,31,32,33,34], soft computing methods [35,36,37,38,39,40,41,42,43,44], and statistical approaches in interaction with fuzzy logic [21,24,45,46,47]. Step 2 is an essential part of the predictive power of the presented model. Fuzzy time series models that rely on imprecise observations have attracted substantial attention in recent years, mainly due to their high applicability to real-life problems.

In fact, a lot of researchers have focused on time series models using imprecise observations. The soft computing techniques employed in this framework are mostly combinations of artificial neural networks, evolutionary algorithms, fuzzy and rough sets. These approaches are widely used for crisp or fuzzy forecasts based on crisp past observations such as electricity load, stock index prices and temperature (for a review of these techniques, we refer to [48,49,50,51,52,53,54,55,56]). In addition, various methods combine techniques of time series and fuzzy regression analysis [57]. For some recent advances in fuzzy regression analysis, see [58,59,60,61,62,63].

The reliability of forecasting methods generally requires exact observations in the sample. But there is often only vague information that is given in terms of imprecise quantities. Moreover, there are various real-world problems related to biological, economic, environmental, medical and sociological data where we face inaccurate instead of accurate data. In many real-life applications, e.g., monthly Co${}_2$ emission, annual sea surface temperature or the water level of a lake, conventional observations are often reported as mean values. In such cases, the data obtained are not sufficient informative since some information contained in the range of the data is neglected. To overcome this shortcoming, one alternative would be to report such kind of data as interval valued (comparable to conventional confidence intervals). However, a potential shortcoming of interval-valued data is the fact that all values within the interval have the same importance. To avoid this issue of interval-valued data, the reported data can alternatively be represented with help of fuzzy numbers [64]. These fuzzy quantities can be modeled via experts opinion, or as simple alternative, they can be constructed via a method proposed by Buckley [65]. In this approach, conventional confidence intervals are employed to construct fuzzy numbers around the conventional mean values.

In addition to the abovementioned methods, there are also fuzzy time series models that rely on fuzzy data, but comparatively few overall. In this regard, Hesamian and Akbari [66] first suggested a fuzzy semi-parametric time series model (FSPTSM) based on fuzzy data, non-fuzzy coefficients, and fuzzy smooth functions. Secondly, Zarei et al. [67] used a specific variant of the FSPTSM [66] for triangular fuzzy data and different distance measures for fuzzy data. And thirdly, Hesamian et al. [68] introduced a forward additive time series model (FATSM) for fuzzy observations.

In this paper we develop a fuzzy nonparametric time series model (FNPTSM) for fuzzy observations that is inspired by nonparametric regression models and kernel smoothing methods [57]. As an initial idea, note that in nonparametric regression analysis, the Nadaraya-Watson estimator [69,70] is fairly common. Now, let us consider the issue of parameter estimation in the nonlinear regression model $x_t=f(x_{t-1},x_{t-2},\dots ,x_{t-p})+{ϵ}_t$ with $f:{\mathbb{R}}^p\to \mathbb{R}$. Based on this general model, a simple nonparametric way of estimating the function f is to employ the kernel-based Nadaraya–Watson estimator

$$\begin{array}{c}\hfill \widehat{f}\left(x_t\right)=\sum _{j=p+1}^{T^{\ast }}{w}^h(t,j)x_j\end{array}$$

(1)

with

$$\begin{array}{c}\hfill w^h(t,j)=\frac{{\sum }_{i=1}^{p}K\left(\frac{x_{t-i}-x_{j-i}}{h}\right)}{{\sum }_{j=p+1}^{T^{\ast }}{\sum }_{i=1}^{p}K\left(\frac{x_{t-i}-x_{j-i}}{h}\right)},\end{array}$$

where K is a kernel function and $h>0$ the bandwidth parameter. Note that the estimator (1) is a weighted average of $x_1,x_2,\dots ,x_T$ using the weights $w^h(t,j)$. As for determining the optimal bandwidth h, the Generalized Cross Validation (GCV) criterion

$$\widehat{h}=arg\underset{h>0}{min}\mathrm{GCV}\left(h\right)=arg\underset{h>0}{min}\frac{1}{T^{\ast }-p}\sum _{t=p+1}^{T^{\ast }}{\left(\frac{x_t-{\sum }_{j=p+1}^{T^{\ast }}{w}^h(t,j)x_j}{1-\frac{\mathrm{tr}\left(W_h\right)}{T^{\ast }-p}}\right)}^2$$

can be utilized, where $\mathrm{tr}\left(W_h\right)$ is the trace of the matrix $W_h=\left[w^h(t,j)\right]$. It is a matter of fact that the estimated values of f are ensured to be within the range of the response variable. This beneficial property is one of the reasons why we apply the Nadaraya–Watson kernel-based estimator for our fuzzy time series model. By utilizing this idea, the proposed FNPTSM provides an estimation procedure of the unknown (nonlinear) relationship between the fuzzy observations in three stages. The advantage of this methodology is that it considerably decreases the complexity in the estimation procedure. While the other fuzzy time series models [66,67,68] are based on estimating the unknown components of the model by unifying the centers and spreads of fuzzy data and their corresponding predicted values, our proposed method provides a smooth estimation procedure according to three separate stages. In the framework of a simulation study and two real-data examples, the efficiency and appropriateness of the FNPTSM is assessed in comparison with previous time series models for fuzzy data by utilizing four approved goodness-of-fit criteria.

The paper is organized as follows. First, we recall some necessary concepts related to fuzzy numbers in Section 2. In Section 3, the three-stage nonparametric kernel-based time series model using fuzzy data is presented. In Section 4, various application examples are given. Concluding remarks are provided in Section 5.

2. Fuzzy Numbers

In this section, we introduce basic definitions of fuzzy numbers that are needed to develop our proposed method.

A fuzzy set $\tilde{A}$ is a mapping on $\mathbb{X}$ that assigns a specific degree of membership $0\le {\mu }_{\tilde{A}}\left(x\right)\le 1$ to each $x\in \mathbb{X}$. In addition, a fuzzy number (FN) $\tilde{A}$ is a convex normalized fuzzy set on the real line $\mathbb{R}$ with an upper semi-continuous membership function of bounded support [71]. In many real applications, vague data a can be reported as $\tilde{A}$: “about a”. Such fuzzy data can often be represented via a special case of FNs, so called $LR$-FNs, which split ${\mu }_{\tilde{A}}$ into two curves: a part on the left and a part on the right of the modal value. So, when considering real-life applications in fuzzy environments, $LR$-FNs play an important role. The membership function of an $LR$-FN${\mu }_{\tilde{A}}\left(x\right)={(a;l_a,r_a)}_{LR}$ can be defined by:

$$\begin{array}{c}\hfill {\mu }_{\tilde{A}}\left(x\right)=\left\{\begin{array}{cc}L\left({\displaystyle \frac{a-x}{l_a}}\right)\hfill & \mathrm{if}x\le a\hfill \\ R\left({\displaystyle \frac{x-a}{r_a}}\right)\hfill & \mathrm{if}x>a\hfill \end{array}\right.\end{array}$$

(2)

In (2), L and R are continuous and strictly decreasing functions from $[0,1]$ to $[0,1]$ satisfying $L\left(0\right)=R\left(0\right)=1$ and $L\left(1\right)=R\left(1\right)=0$. In addition, $a\in \mathbb{R}$ represents the modal value, while $l_a>0$ and $r_a>0$ are the left spreads and right spreads of $\tilde{A}$, respectively. The set of all $LR$-FNs is represented by ${\mathcal{F}}_{LR}\left(\mathbb{R}\right)$. A special case of an $LR$-FN is the so-called triangular fuzzy number (TFN), whose membership function has the following form:

$${\mu }_{\tilde{A}}\left(x\right)=\left\{\begin{array}{cc}\frac{x-(a-l_a)}{l_a}\hfill & a-l_a\le x\le a\hfill \\ \frac{a+r_a-x}{r_a}\hfill & a<x\le a+r_a\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

There are various operations that can be defined between two $LR$-FNs, i.e., between $\tilde{A}={(a;l_a,r_a)}_{LR}$ and $\tilde{B}={(b;l_b,r_b)}_{LR}$. For instance, as we need both operations in this paper, we define Addition and Scalar multiplication of $\tilde{A}$ and $\tilde{B}$ in the following [72]:

• Addition: $\tilde{A}\oplus \tilde{B}={(a+b;l_a+l_b,r_a+r_b)}_{LR}$
• Scalar multiplication:

  $$\lambda \otimes \tilde{A}=\left\{\begin{array}{cc}{(\lambda a;\lambda l_a,\lambda r_a)}_{LR}\hfill & \mathrm{if}\lambda >0\hfill \\ {(\lambda a;-\lambda r_a,-\lambda l_a)}_{RL}\hfill & \mathrm{if}\lambda <0\hfill \end{array}\right.$$

Moreover, there are numerous concepts used to define distances between two $LR$-FNs $\tilde{A}={(a;l_a,r_a)}_{LR}$ and $\tilde{B}={(b;l_b,r_b)}_{LR}$ [71]. Here, we utilize the squared error distance measure D for performance evaluation of the FNPTSM in comparison with other models. It is defined as

$$D(\tilde{A},\tilde{B})={(({(a-b)}^2+c_1{(l_a-l_b)}^2+c_2{(r_a-r_b)}^2)/3)}^{0.5}$$

with $c_1={\int }_0^1{L}^{-1}\left(\alpha \right)d\alpha$ and $c_2={\int }_0^1{R}^{-1}\left(\alpha \right)d\alpha$ [73].

3. Nonparametric Kernel-Based Time Series Model for Fuzzy Data

In this section, the FNPTSM is developed along with the suggested parameter estimation method.

3.1. The Model

First, we recall the definition of fuzzy time series data.

Definition 1. 

Let ${\tilde{\mathit{x}}}_{\mathit{T}}=\{{\tilde{x}}_1,{\tilde{x}}_2,\dots ,{\tilde{x}}_T\}$ be a set of FNs of size T. Then, ${\tilde{\mathit{x}}}_{\mathit{T}}$ is called fuzzy time series data if $\{{\tilde{x}}_1,{\tilde{x}}_2,\dots ,{\tilde{x}}_T\}$ is the vague concept of ordinary time series data $\{x_1,x_2,\dots ,x_T\}$ [68,74].

As discussed in the Introduction, there are many situations where it is preferable to report exact data x by an FN$\tilde{x}$ as “about x”. Then, $\tilde{x}$ is the respective vague concept of x.

Definition 2. 

Let ${\tilde{\mathit{x}}}_{\mathit{T}}=\{{\tilde{x}}_1,{\tilde{x}}_2,\dots ,{\tilde{x}}_T\}$ be fuzzy time series data. The FNPTSM for fuzzy time series data ${\tilde{\mathit{x}}}_{\mathit{T}}$ is then defined by

$${\tilde{x}}_t=\tilde{f}({\tilde{x}}_{t-1},{\tilde{x}}_{t-2},\dots ,{\tilde{x}}_{t-p})\oplus {\widetilde{ϵ}}_t,$$

(3)

where

• ${\tilde{x}}_t={(x_t;l_{x_t},r_{x_t})}_{LR}$,
• $\tilde{f}({\tilde{x}}_{t-1},{\tilde{x}}_{t-2},\dots ,{\tilde{x}}_{t-p})=(f(x_{t-1},x_{t-2},\dots ,x_{t-p});$
  $l_{f(l_{x_{t-1}},l_{x_{t-2}},\dots ,l_{x_{t-p}})},r_{f(r_{x_{t-1}},r_{x_{t-2}},\dots ,r_{x_{t-p}})}{)}_{LR}$,
• ${\widetilde{ϵ}}_t={({ϵ}_t;l_{{ϵ}_t},r_{{ϵ}_t})}_{LR}$’s are fuzzy errors, where ${ϵ}_t\in \mathbb{R}$ and $l_{{ϵ}_t},r_{{ϵ}_t}\in {\mathbb{R}}^{+}$.

Remark 1. 

Note that (3) provides an FN in the form ${\tilde{x}}_t^{\ast }=$$(x_t^{\ast };l_{x_t^{\ast }},$$r_{x_t^{\ast }}{)}_{LR}$ with $x_t^{\ast }=$$f(x_{t-1},x_{t-2},$$\dots ,x_{t-p})+{ϵ}_t$, $l_{x_t^{\ast }}=$$l_{f(l_{x_{t-1}},\dots ,l_{x_{t-p}})+l_{{ϵ}_t}}$ and $r_{x_t^{\ast }}=r_{f(r_{x_{t-1}},\dots ,r_{x_{t-p}})+r_{{ϵ}_t}}$ with $t=1,2,\dots ,T$. According to Definition 1, as $\{x_1,x_2,\dots ,x_T\}$ is ordinary time series data, ${\tilde{\mathit{x}}}_{\mathit{T}}^{\ast }=\{{\tilde{x}}_1^{\ast },{\tilde{x}}_2^{\ast },\dots ,{\tilde{x}}_T^{\ast }\}$ is also a vague concept of ordinary time series data $\{x_1^{\ast },x_2^{\ast },\dots ,x_T^{\ast }\}$. Thus, the proposed fuzzy time series model (3) generates new fuzzy time series data.

3.2. Three-Stage Estimation Method for the Nonlinear Fuzzy Smooth Function

Below, we suggest a three-stage method to estimate the unknown fuzzy smooth function $\tilde{f}$ in (3). For this purpose, the fuzzy predictions are obtained based on a within-sample forecast ${\mathit{x}}_{T^{\ast }}={\left(x_1,x_2,\dots ,x_{T^{\ast }}\right)}^{\top }$ with $T^{\ast }<T$. From (3), one can get three ordinary nonlinear time series models as (1) $x_t=f(x_{t-1},x_{t-2},\dots ,x_{t-p})+{ϵ}_t$, (2) $l_{x_t}=l_{f(l_{x_{t-1}},\dots ,l_{x_{t-p}})+l_{{ϵ}_t}}$, and (3) $r_{x_t}=r_{f(r_{x_{t-1}},\dots ,r_{x_{t-p}})+r_{{ϵ}_t}}$ for $t=1,2,\dots ,T^{\ast }$. Therefore, to estimate the fuzzy smooth function at $\tilde{\mathit{x}}=$${(\mathit{x};l_{\mathit{x}},r_{\mathit{x}})}_T$ with $\mathit{x}=$${(x_1,x_2,\dots ,x_p)}^{\top }$, $l_{\mathit{x}}=$$(l_{x_1},$$l_{x_2},\dots ,$$l_{x_p}{)}^{\top }$ and $r_{\mathit{x}}=$$(r_{x_1},$$r_{x_2},\dots ,$$r_{x_p}{)}^{\top }$, we follow the three-stage procedure below:

• Stage (1): Consider the nonlinear regression model $l_{x_t}=l_{f(x_{t-1},x_{t-2},\dots ,x_{t-p})}+l_{{ϵ}_t}$. Based on the time series data $l_{{\mathit{x}}_t}={(l_{x_{t-1}},\dots ,l_{x_{t-p}})}^{\top }$, we employ the weighted Nadaraya–Watson estimator to estimate $l_f$ for a within-sample forecast $T^{\ast }\le T$ at $l_{\mathit{x}}={(l_{x_1},\dots ,l_{x_p})}^{\top })$ as

  $$l_{\widehat{f}\left(l_{{\mathit{x}}_t}\right)}=\sum _{j=p+1}^{T^{\ast }}{w}^{h_l}(t,j)l_{x_j},$$

  where

  $$\begin{array}{c}\hfill w^{h_l}(t,j)=\frac{{\sum }_{i=1}^{p}K\left(\frac{l_{x_{t-i}}-l_{x_{j-i}}}{h_l}\right)}{{\sum }_{j=p+1}^{T^{\ast }}{\sum }_{i=1}^{p}K\left(\frac{l_{x_{t-i}}-l_{x_{j-i}}}{h_l}\right)}\end{array}$$

  (4)

  with kernel function $K(.)$ and bandwidth parameter $h_l>0$. The optimal value of $h_l$ can be estimated by implementing the GCV criterion,

  $${\widehat{h}}_l=arg\underset{h_l>0}{min}\mathrm{GCV}\left(h\right)=arg\underset{h_l>0}{min}\frac{1}{T^{\ast }-p}\sum _{t=p+1}^{T^{\ast }}{\left(\frac{l_{x_t}-{\sum }_{j=p+1}^{T^{\ast }}{w}^{h_l}(t,j)l_{x_j}}{1-\frac{\mathrm{tr}\left(W_{h_l}\right)}{T^{\ast }-p}}\right)}^2,$$

  (5)

  where $\mathrm{tr}\left(W_{h_l}\right)$ is the trace of the matrix $W_{h_l}=\left[w^{h_l}(t,j)\right]$ with $w^{h_l}(t,j)$ as defined in (4).
• Stage (2): Consider the nonlinear regression model $r_{x_t}=r_{f(x_{t-1},x_{t-2},\dots ,x_{t-p})}+r_{{ϵ}_t}$. Based on the within-sample time series forecast data $r_{{\mathit{x}}_t}={(r_{x_{t-1}},\dots ,r_{x_{t-p}})}^{\top }$, $t=1,2,\dots ,T^{\ast }$, the weighted Nadaraya–Watson estimation of $r_f$ at $r_{\mathit{x}}={(r_{x_1},\dots ,r_{x_p})}^{\top })$ can be established via

  $$r_{\widehat{f}\left(r_{{\mathit{x}}_t}\right)}=\sum _{j=p+1}^{T^{\ast }}{w}^{h_r}(t,j)r_{x_j},$$

  where

  $$\begin{array}{c}\hfill w^{h_r}(t,j)=\frac{{\sum }_{i=1}^{p}K\left(\frac{r_{x_{t-i}}-r_{x_{j-i}}}{h_r}\right)}{{\sum }_{j=p+1}^{T^{\ast }}{\sum }_{i=1}^{p}K\left(\frac{r_{x_{t-i}}-r_{x_{j-i}}}{h_r}\right)}\end{array}$$

  (6)

  and $h_r>0$ is a bandwidth parameter. The optimal value of $h_r$ can be estimated using the GCV criterion,

  $${\widehat{h}}_r=arg\underset{h_r>0}{min}\mathrm{GCV}\left(h\right)=arg\underset{h_r>0}{min}\frac{1}{T^{\ast }-p}\sum _{t=p+1}^{T^{\ast }}{\left(\frac{r_{x_t}-{\sum }_{j=p+1}^{T^{\ast }}{w}^{h_r}(t,j)r_{x_j}}{1-\frac{\mathrm{tr}\left(W_{h_r}\right)}{T^{\ast }-p}}\right)}^2,$$

  (7)

  where $\mathrm{tr}\left(W_{h_r}\right)$ is the trace of the matrix $W_{h_r}=\left[w^{h_r}(t,j)\right]$ with $w^{h_r}(t,j)$ as defined in (6).
• Stage (3): Consider the nonlinear regression model $x_t=$$f(x_{t-1},x_{t-2},\dots ,x_{t-p})+{ϵ}_t$. Based on the within-sample time series forecast data $({\mathit{x}}_t=$$(x_{t-1},$$x_{t-2},\dots ,$$x_{t-p}{)}^{\top })$, $t=1,2,\dots ,T^{\ast }$, a nonparametric estimator f can be achieved as

  $$\widehat{f}\left({\mathit{x}}_t\right)=\sum _{j=p+1}^{T^{\ast }}{w}^h(t,j)x_j,$$

  where

  $$\begin{array}{c}\hfill w^h(t,j)=\frac{{\sum }_{i=1}^{p}K\left(\frac{x_{t-i}-x_{j-i}}{h}\right)}{{\sum }_{j=p+1}^{T^{\ast }}{\sum }_{i=1}^{p}K\left(\frac{x_{t-i}-x_{j-i}}{h}\right)}\end{array}$$

  (8)

  and bandwidth parameter $h>0$. Similar to the previous stages, the optimal value of h is estimated with the help of the GCV criterion,

  $$\widehat{h}=arg\underset{h>0}{min}\mathrm{GCV}\left(h\right)=arg\underset{h>0}{min}\frac{1}{T^{\ast }-p}\sum _{t=p+1}^{T^{\ast }}{\left(\frac{x_t-{\sum }_{j=p+1}^{T^{\ast }}{w}^h(t,j)x_j}{1-\frac{\mathrm{tr}\left(W_h\right)}{T^{\ast }-p}}\right)}^2,$$

  (9)

  where $\mathrm{tr}\left(W_h\right)$ is the trace of the matrix $W_h=\left[w^h(t,j)\right]$ with $w^h(t,j)$, as defined in (8).

Therefore, the forecast ${\tilde{x}}_{T^{\ast }+k}$ with time lag $k\in \mathbb{N}$ can be achieved by an $LR$-FN via ${\tilde{\widehat{x}}}_{T^{\ast }+k}=$$({\widehat{x}}_{T^{\ast }+k};$$l_{{\widehat{x}}_{T^{\ast }+k}},$$r_{{\widehat{x}}_{T^{\ast }+k}}{)}_{LR}$ with

$$\begin{array}{cc}\hfill {\widehat{x}}_{T^{\ast }+k}& =\sum _{j=p+1}^{T^{\ast }+k-1}\frac{{\sum }_{i=1}^{p}K\left(\frac{x_{t-i}-x_{j-i}}{\widehat{h}}\right)}{{\sum }_{j=p+1}^{T^{\ast }+k-1}{\sum }_{i=1}^{p}K\left(\frac{x_{t-i}-x_{j-i}}{\widehat{h}}\right)}·x_j,\hfill \\ \hfill l_{{\widehat{x}}_{T^{\ast }+k}}& =\sum _{j=p+1}^{T^{\ast }+k-1}\frac{{\sum }_{i=1}^{p}K\left(\frac{l_{x_{t-i}}-l_{x_{j-i}}}{{\widehat{h}}_l}\right)}{{\sum }_{j=p+1}^{T^{\ast }+k-1}{\sum }_{i=1}^{p}K\left(\frac{l_{x_{t-i}}-l_{x_{j-i}}}{{\widehat{h}}_l}\right)}·l_{x_j},\hfill \\ \hfill r_{{\widehat{x}}_{T^{\ast }+k}}& =\sum _{j=p+1}^{T^{\ast }+k-1}\frac{{\sum }_{i=1}^{p}K\left(\frac{r_{x_{t-i}}-r_{x_{j-i}}}{{\widehat{h}}_r}\right)}{{\sum }_{j=p+1}^{T^{\ast }+k}{\sum }_{i=1}^{p}K\left(\frac{r_{x_{t-i}}-r_{x_{j-i}}}{{\widehat{h}}_r}\right)}·r_{x_j}.\hfill \end{array}$$

According to Stages (2) and (3), it can be seen that the spreads of the fuzzy prediction ${\tilde{x}}_{T^{\ast }+k}$ are always non-negative.

Remark 2. 

Since the proposed time series model relies on fuzzy data, let us recall the previous time series models based on fuzzy data [66,67,68]. First, Hesamian and Akbari [66] proposed a fuzzy semi-parametric autoregressive integrated moving average (ARIMA) model as follows:

$${\tilde{x}}_i=\underset{l=1}{\overset{p}{⨁}}({\theta }_l\otimes {\tilde{x}}_{i-l}\oplus \tilde{f}\left(t_i\right)\oplus {\widetilde{ϵ}}_i),i=p+1,\dots ,T.$$

The parameters of the model are estimated by employing a hybrid method including a nonparametric kernel-based method and least absolute deviations. For a second time series model based on fuzzy data, Zarei et al. [67] applied the method [66] to estimate the model parameters and the fuzzy smooth function based on a specific distance, kernel and triangular fuzzy numbers. Finally, Hesamian et al. [68] proposed the fuzzy nonlinear time series model

$${\tilde{x}}_t=\tilde{f}({\tilde{x}}_{t-1},{\tilde{x}}_{t-2},\dots ,{\tilde{x}}_{t-p})\oplus {\widetilde{ϵ}}_t,t=1,2,\dots ,T,$$

where

$$\tilde{f}({\tilde{x}}_{t-1},{\tilde{x}}_{t-2},\dots ,{\tilde{x}}_{t-p})=\underset{l=1}{\overset{p}{⨁}}{f}_l\left({\tilde{x}}_{t-l}\right).$$

As for the estimation of the unknown fuzzy smooth functions ${\tilde{f}}_l$, they applied a forward additive nonparametric technique.

Remark 3. 

We have extended some common performance measures used to compare the predictive accuracy of different time series models that we implement in Section 4. For this purpose, a time series model is first estimated based on a within-sample fuzzy time series dataset of size $T^{\ast }<T$ and then the performance of the model is evaluated via the remaining fuzzy time series dataset of size $T-T^{\ast }$.

• Mean Forecast Error:

  $$\mathit{MFE}=\frac{{\sum }_{t=T^{\ast }+1}^T{D}^2({\tilde{\widehat{x}}}_t,{\tilde{x}}_t)}{T-T^{\ast }}$$
• Mean Absolute Scaled Error:

  $$\mathit{MASE}=\frac{{\sum }_{t=T^{\ast }+1}^T{q}_t}{T-T^{\ast }}$$
  with

  $$q_t=\frac{D({\tilde{\widehat{x}}}_t,{\tilde{x}}_t)}{\frac{1}{T-T^{\ast }}{\sum }_{t=T^{\ast }+1}^T{D}^2({\tilde{x}}_t,{\tilde{x}}_{t-1})}$$
• Basis of the Index of Agreement:

  $$\mathit{BIA}=1-{\displaystyle \frac{{\sum }_{t=T^{\ast }+1}^T{D}^2({\tilde{x}}_t,{\tilde{\widehat{x}}}_t)}{{\sum }_{t=T^{\ast }+1}^T{(D({\tilde{x}}_t,\tilde{\overline{x}})+D(\tilde{\overline{x}},{\tilde{\widehat{x}}}_t))}^2}}$$
  with

  $$\tilde{\overline{x}}=\frac{{\sum }_{t=T^{\ast }+1}^T{\tilde{x}}_t}{T-T^{\ast }}$$
• Mean Similarity Measure:

  $$\mathit{MSM}=\frac{1}{T-T^{\ast }}\sum _{t=T^{\ast }+1}^T\frac{\int min\{{\tilde{\widehat{x}}}_t\left(x\right),{\tilde{x}}_t\left(x\right)\}dx}{\int max\{{\tilde{\widehat{x}}}_t\left(x\right),{\tilde{x}}_t\left(x\right)\}dx}$$

Let A and B be two fuzzy time series models. As $\mathit{MSM}:{\mathcal{F}}_{LR}\left(\mathbb{R}\right)\times {\mathcal{F}}_{LR}\left(\mathbb{R}\right)\to [0,1]$ is a similarity measure, values of MSM above 0.5 show a good degree of similarity between the fuzzy responses and their fuzzy predictions. If we observe ${\mathit{MSM}}_B<{\mathit{MSM}}_A$, then model A outperforms model B. Further, if ${\mathit{MFE}}_A<{\mathit{MFE}}_B$, ${\mathit{MASE}}_A<{\mathit{MASE}}_B$ or ${\mathit{BIA}}_A<{\mathit{BIA}}_B$, then model A acts better in terms of prediction accuracy compared to model B.

Remark 4. 

While the proposed estimation procedure does not depend on the shape functions L and R corresponding to fuzzy data, the performance measures MFE, MASE and MSM depend on these shape functions. Therefore, the selected type of the shape functions L and R may affect the prediction criteria. For instance, assume that the data have reported by ${\tilde{x}}_t={(x_t,l_{x_t},r_{x_t})}_{LR}$ with $L\left(x\right)=1-x$ and $R\left(x\right)=\sqrt{1-x}$. That is, $c_1=\frac{1}{2}$ and $c_2=\frac{2}{3}$. Therefore, the distance between ${\tilde{x}}_t$ and its prediction is $D^2({\tilde{x}}_t,{\tilde{\widehat{x}}}_t)={(x_t-{\widehat{x}}_t)}^2+\frac{1}{2}{(l_{x_t}-l_{{\widehat{x}}_t})}^2+\frac{2}{3}{(r_{x_t}-r_{{\widehat{x}}_t})}^2$. This implies that the MFE criterion is more sensitive to right spreads than to left spreads in this case. Considering $L\left(x\right)=R\left(x\right)=1-x$, it can be seen that $D^2({\tilde{x}}_t,{\tilde{\widehat{x}}}_t)$ would be equally dependent from the left and right spreads. However, when we compare the performance of fuzzy time series models, it is reasonable that the shape functions L and R are assumed to be the same for all the considered models. Thus, following this approach, the performance criteria are not sensitive to the selection of L and R since $c_1$ and $c_2$ remain fixed for each model.

3.3. Selection of Autoregressive Order and Optimal Bandwidths

When implementing the FNPTSM (3), it is necessary to select the optimal bandwidths h, $h_l$ and $h_r$, to choose the kernel function and to determine the autoregressive order p. The procedure used to select the autoregressive order and the optimal bandwidths is proposed as follows:

(1)

Let $p=1$.

(2)

(2.1) Compute ${\widehat{h}}_l^p$ based on (5).

(2.2) Compute ${\widehat{h}}_r^p$ based on (7).

(2.3) Compute ${\widehat{h}}^p$ based on (9).

(3)

Let $p=p+1$ and return to (2) until

$$\widehat{p}=arg\underset{p}{min}{\mathrm{RMSE}}_p,$$

where

$${\mathrm{RMSE}}_p=\sqrt{\frac{{\sum }_{i=p+1}^{T^{\ast }}{D}^2({\tilde{\widehat{x}}}_i,{\tilde{x}}_i)}{T^{\ast }-p}}.$$

Then, $\widehat{p}$, ${\widehat{h}}_p$, ${\widehat{h}}_l^p$ and ${\widehat{h}}_r^p$ are the optimal values.

4. Numerical Examples

In this section, the effectiveness of the FNPTSM is investigated considering a simulation study and application examples that rely on fuzzy data. Recall that there are three other time series models that are based on fuzzy data (see Remark 2), i.e., the models introduced by Hesamian and Akbari [66], Zarei et al. [67] and Hesamian et al. [68]. However, as the method of Zarei et al. [67] is based on Hesamian and Akbari’s method [66] (with a different distance measure), we omit this technique in the comparisons below. Thus, we compare our proposed method with the models suggested by Hesamian and Akbari (FSPTSM) [66] and Hesamian et al. (FATSM) [68] via three different kernel functions (Gaussian, Epanechnikov, and triweight).

Example 1. 

In this example, 10 fuzzy datasets, each of size 300, are generated by the following FNPTSM:

$${\tilde{x}}_t=f({\tilde{x}}_{t-1},{\tilde{x}}_{t-2},{\tilde{x}}_{t-3})\oplus {\widetilde{ϵ}}_t,t=4,5,\dots ,300,$$

where

1.

$$f({\tilde{x}}_1,{\tilde{x}}_2,{\tilde{x}}_3)={\left(x_1-cos\left(x_2\right)-exp\left(\frac{x_3}{1+|x_3|}\right);{cos}^2\left(0.9\prod _{j=1}^3{l}_{x_j}\right),exp\left(0.002\prod _{j=1}^3{r}_{x_j}\right)\right)}_{LR}$$

2.

${\tilde{x}}_j={(x_j;l_{x_j},r_{x_j})}_{LR}$, $j=1,2,3$ are the initial values with $x_j\sim N(0,1)$, and $l_{x_j}$ and $r_{x_j}$ are random variables following $U(0,0.2)$ and $U(0,0.9)$, respectively,

3.

${\widetilde{ϵ}}_t={({ϵ}_t;l_{{ϵ}_t},r_{{ϵ}_t})}_{LR}$ with ${ϵ}_t\sim N(0,4)$, $l_{{ϵ}_t}$ and $r_{{ϵ}_t}$ are random variables following $U(0,0.4)$ and $U(0,0.5)$, respectively, and

4.

$L\left(x\right)=1-x^2$ and $R\left(x\right)=1-x$.

The kernels Gaussian, Epanechnikov, and triweight are applied to predict ${\tilde{x}}_t$. Based on the 10 sample fuzzy datasets (each of size 300), the mean values of the goodness-of-fit measures and their corresponding bandwidth mean values are summarized in

Table 1. The mean performance measures of the FNPTSM, FSPTSM and FATSM corresponding to some specific kernels in Example 1.

Method	Kernel	Results	Goodness-of-Fit Criteria
FNPTSM	Gaussian		$\overline{\mathrm{MFE}}=1.0452$
		$\overline{\widehat{h}}=0.45$	$\overline{\mathrm{MASE}}=1.6089$
		$\overline{\widehat{h^l}}=0.04$	$\overline{\mathrm{BIA}}=0.9996$
		$\overline{\widehat{h^r}}=0.22$	$\overline{\mathrm{MSM}}=0.4167$
	Epanechnikov		$\overline{\mathrm{MFE}}=1.1728$
		$\overline{\widehat{h}}=0.66$	$\overline{\mathrm{MASE}}=1.6478$
		$\overline{\widehat{h^l}}=0.05$	$\overline{\mathrm{BIA}}=0.9992$
		$\overline{\widehat{h^r}}=0.39$	$\overline{\mathrm{MSM}}=0.3953$
	triweight		$\overline{\mathrm{MFE}}=1.2482$
		$\overline{\widehat{h}}=1.89$	$\overline{\mathrm{MASE}}=1.6339$
		$\overline{\widehat{h^l}}=0.08$	$\overline{\mathrm{BIA}}=0.9991$
		$\overline{\widehat{h^r}}=0.54$	$\overline{\mathrm{MSM}}=0.3721$
FSPTSM	Gaussian	${\overline{h}}_{\mathrm{opt}}=0.07$	$\overline{\mathrm{MFE}}=8.9728$
		${\widehat{\theta }}_1=0.5575$	$\overline{\mathrm{MASE}}=4.2652$
		${\widehat{\theta }}_2=-0.0956$	$\overline{\mathrm{BIA}}=0.9536$
		${\widehat{\theta }}_3=-0.1247$	$\overline{\mathrm{MSM}}=0.2207$
	Epanechnikov	${\overline{h}}_{\mathrm{opt}}=0.02$	$\overline{\mathrm{MFE}}=5.2530$
		${\widehat{\theta }}_1=-0.6424$	$\overline{\mathrm{MASE}}=4.5383$
		${\widehat{\theta }}_2=-0.5168$	$\overline{\mathrm{BIA}}=0.9603$
		${\widehat{\theta }}_3=-0.4258$	$\overline{\mathrm{MSM}}=0.2920$
	triweight	${\overline{h}}_{\mathrm{opt}}=0.13$	$\overline{\mathrm{MFE}}=11.5231$
		${\widehat{\theta }}_1=0.3757$	$\overline{\mathrm{MASE}}=3.8643$
		${\widehat{\theta }}_2=-0.1576$	$\overline{\mathrm{BIA}}=0.9738$
		${\widehat{\theta }}_3=-0.2568$	$\overline{\mathrm{MSM}}=0.2133$
FATSM	Gaussian		$\overline{\mathrm{MFE}}=1.392$
		${\overline{\widehat{h}}}_1=1.8$	$\overline{\mathrm{MASE}}=18.028$
		${\overline{\widehat{h}}}_2=0.7$	$\overline{\mathrm{BIA}}=0.972$
		${\overline{\widehat{h}}}_3=0.03$	$\overline{\mathrm{MSM}}=0.321$
	Epanechnikov		$\overline{\mathrm{MFE}}=1.353$
		${\overline{\widehat{h}}}_1=1.75$	$\overline{\mathrm{MASE}}=17.548$
		${\overline{\widehat{h}}}_2=1.20$	$\overline{\mathrm{BIA}}=0.976$
		${\overline{\widehat{h}}}_3=0.6$	$\overline{\mathrm{MSM}}=0.339$
	triweight		$\overline{\mathrm{MFE}}=1.348$
		${\overline{\widehat{h}}}_1=2$	$\overline{\mathrm{MASE}}=16.459$
		${\overline{\widehat{h}}}_2=1.5$	$\overline{\mathrm{BIA}}=0.979$
		${\overline{\widehat{h}}}_3=0.5$	$\overline{\mathrm{MSM}}=0.349$

. Consulting the results for the FNPTSM, it is evident that the best results among various kernels are obtained via the Gaussian kernel (lowest values of $\overline{\mathit{MFE}}$, $\overline{\mathit{MASE}}$ and largest values of $\overline{\mathit{BIA}}$, $\overline{\mathit{MSM}}$). In addition, the results of the FSPTSM and FATSM can also be found in Table 1. Comparing these results with the results of the FNPTSM, it is obvious that the FNPTSM provides more accurate predictions compared to both other methods for all three kernels, as all the considered goodness-of-fit measures show better results for the FNPTSM. That is, we observe the lowest values of $\overline{\mathit{MFE}}$, $\overline{\mathit{MASE}}$ and the largest values of $\overline{\mathit{BIA}}$, $\overline{\mathit{MSM}}$ for the FNPTSM.

Example 2. 

Three models, the FNPTSM, FSPTSM and FATSM, are implemented to analyze the dataset in

Table 2. Fuzzy time series data in Example 2.

t	${\tilde{\mathit{x}}}_{\mathit{t}}$	t	${\tilde{\mathit{x}}}_{\mathit{t}}$
1	${(1.7337;0.8051)}_T$	15	${(2.9145;1.1507)}_T$
2	${(2.3302;0.9228)}_T$	16	${(2.6085;1.1335)}_T$
3	${(1.3199;0.7742)}_T$	17	${(3.0432;0.4489)}_T$
4	${(5.0507;0.8948)}_T$	18	${(6.8010;0.9588)}_T$
5	${(1.4206;1.0540)}_T$	19	${(4.9351;0.8115)}_T$
6	${(4.0273;0.9331)}_T$	20	${(3.5672;0.6054)}_T$
7	${(2.8624;1.0480)}_T$	21	${(3.8828;1.1579)}_T$
8	${(4.7107;1.0647)}_T$	22	${(0.5183;0.9652)}_T$
9	${(4.1098;1.1028)}_T$	23	${(3.6846;0.8175)}_T$
10	${(4.4843;1.0881)}_T$	24	${(3.5117;0.4248)}_T$
11	${(1.3249;1.0064)}_T$	25	${(2.8294;0.7956)}_T$
12	${(3.2249;0.5503)}_T$	26	${(2.3836;1.0101)}_T$
13	${(3.2916;0.5049)}_T$	27	${(3.9454;0.7558)}_T$
14	${(3.5508;0.9268)}_T$	28	${(3.6012;1.0266)}_T$

taken from [67].

Eighty percent of the data were used for parameter estimation and the rest were applied to fit the model. The goodness-of-fit values are given in

Table 3. The performance measures of the FNPTSM, FSPTSM and FATSM corresponding to some specific kernels in Example 2.

Method	Kernel	Results	Goodness-of-Fit Criteria
FNPTSM	Gaussian	$\widehat{p}=2$	$\mathrm{MFE}=0.1127$
		$\widehat{h}=0.64$	$\mathrm{MASE}=0.4030$
		$\widehat{h^l}=0.22$	$\mathrm{BIA}=0.9629$
		$\widehat{h^r}=0.22$	$\mathrm{MSM}=0.3791$
	Epanechnikov	$\widehat{p}=2$	$\mathrm{MFE}=0.1133$
		$\widehat{h}=1.35$	$\mathrm{MASE}=0.4016$
		$\widehat{h^l}=0.09$	$\mathrm{BIA}=0.9639$
		$\widehat{h^r}=0.09$	$\mathrm{MSM}=0.3836$
	triweight	$\widehat{p}=3$	$\mathrm{MFE}=0.1087$
		$\widehat{h}=1.44$	$\mathrm{MASE}=0.3722$
		$\widehat{h^l}=0.65$	$\mathrm{BIA}=0.9626$
		$\widehat{h^r}=0.65$	$\mathrm{MSM}=0.4439$
FSPTSM	Gaussian	$\widehat{h}=0.13,\widehat{p}=3$	$\mathrm{MFE}=3.6124$
		${\widehat{\theta }}_1=-0.2305$	$\mathrm{MASE}=2.1452$
		${\widehat{\theta }}_2=-0.0613$	$\mathrm{BIA}=0.9666$
		${\widehat{\theta }}_3=-0.4164$	$\mathrm{MSM}=0.0061$
	Epanechnikov	$\widehat{h}=0.16$	$\mathrm{MFE}=1.5749$
		$\widehat{p}=1$	$\mathrm{MASE}=1.3135$
		${\widehat{\theta }}_1=-0.3192$	$\mathrm{BIA}=0.9438$
			$\mathrm{MSM}=0.1080$
	triweight	$\widehat{h}=0.23$	$\mathrm{MFE}=1.6252$
		$\widehat{p}=1$	$\mathrm{MASE}=1.3371$
		${\widehat{\theta }}_1=-0.3437$	$\mathrm{BIA}=0.9443$
			$\mathrm{MSM}=0.1205$
FATSM	Gaussian	$\widehat{p}=2$	$\mathrm{MFE}=0.327$
		${\widehat{h}}_1=0.2$	$\mathrm{MASE}=0.943$
		${\widehat{h}}_2=0.3$	$\mathrm{BIA}=0.954$
			$\mathrm{MSM}=0.341$
	Epanechnikov	$\widehat{p}=2$	$\mathrm{MFE}=0.386$
		${\widehat{h}}_1=0.3$	$\mathrm{MASE}=1.302$
		${\widehat{h}}_2=1.5$	$\mathrm{BIA}=0.948$
			$\mathrm{MSM}=0.3208$
	triweight	$\widehat{p}=3$	$\mathrm{MFE}=0.254$
		${\widehat{h}}_1=0.5$	$\mathrm{MASE}=0.733$
		${\widehat{h}}_2=0.1$	$\mathrm{BIA}=0.958$
		${\widehat{h}}_3=0.05$	$\mathrm{MSM}=0.358$

for the three kernels. The best results among various kernels are obtained by employing the triweight kernel (lowest values of MFE, MASE, and largest value of MSM). The results of the FSPTSM and FATSM are also given in Table 3. As for the FSPTSM , the best results are obtained based on the triweight kernel with $\mathit{MFE}=1.6252$, $\mathit{MASE}=1.3371$, $\mathit{MSM}=0.1205$ and $\mathit{BIA}=0.9443$. The best results of FATSM are also obtained based on the triweight kernel with $\mathit{MFE}=0.254$, $\mathit{MASE}=0.733$, $\mathit{MSM}=0.358$ and $\mathit{BIA}=0.958$. However, all goodness-of-fit measures related to the FNPTSM show a better performance compared to both the FSPTSM and FATSM, i.e., the lowest values of MFE, MASE and the largest values of BIA, MSM are observed for the FNPTSM. The results show that the newly presented FNPTSM is more efficient than the FSPTSM and FATSM for the data in Table 3. The plot of the fuzzy data and corresponding estimates based on the triweight kernel is given in Figure 1 for all methods (FNPTSM , FSPTSM, FATSM).

Example 3. 

In this example, we employ the FNPTSM and both the FSPTSM and FATSM to predict the global land–ocean temperature [75]. For this purpose, we use the global land–ocean temperature from January 2000 to December 2020, as shown in Figure 2. The data are reported as average values for each month. Therefore, the data can also be interpreted as “mean of each month" and appropriately modeled via triangular fuzzy numbers. Inspired by Buckley [65], this dataset can be used to evaluate the global land-ocean temperature with the help of a TFN ${\tilde{x}}_t={({\overline{x}}_t;Z_{0.005}{s}_t/\sqrt{n_t},0.15Z_{0.025}{s}_t/\sqrt{n_t})}_T$, where $n_t$, ${\overline{x}}_t$, and $s_t$ denote the number of days, mean, and standard deviation of the global land–ocean temperature in the $t^{\mathrm{th}}$ month, respectively, and $Z_{\alpha }$ is the α-quantile of the standard normal distribution. However, since we do not have daily values of the global land–ocean temperature, we model the monthly global land–ocean temperature for a month t via ${\tilde{x}}_t$$=(x_t;$$0.2x_t,$$0.15x_t{)}_T$.

In this example, 200 observations were used to estimate the parameters. A further 52 observations were used to fit the model. The goodness-of-fit values that correspond to the FNPTSM, FSPTSM and FATSM are given in

Table 4. The performance measures of the FNPTSM, FSPTSM and FATSM corresponding to some specific kernels in Example 3.

Method	Kernel	Results	Goodness-of-Fit Criteria
FNPTSM	Gaussian	$\widehat{p}=2$	$\mathrm{MFE}=0.0061$
		$\widehat{h}=0.030$	$\mathrm{MASE}=14.6082$
		$\widehat{h^l}=0.006$	$\mathrm{BIA}=0.9783$
		$\widehat{h^r}=0.005$	$\mathrm{MSM}=0.3781$
	Epanechnikov	$\widehat{p}=2$	$\mathrm{MFE}=0.0061$
		$\widehat{h}=0.051$	$\mathrm{MASE}=14.6360$
		$\widehat{h^l}=0.013$	$\mathrm{BIA}=0.9782$
		$\widehat{h^r}=0.007$	$\mathrm{MSM}=0.3787$
	triweight	$\widehat{p}=2$	$\mathrm{MFE}=0.0067$
		$\widehat{h}=0.032$	$\mathrm{MASE}=14.8885$
		$\widehat{h^l}=0.011$	$\mathrm{BIA}=0.9761$
		$\widehat{h^r}=0.010$	$\mathrm{MSM}=0.3921$
FSPTSM	Gaussian	$\widehat{h}=0.05,\widehat{p}=3$	$\mathrm{MFE}=0.0153$
		${\widehat{\theta }}_1=0.3912$	$\mathrm{MASE}=20.4387$
		${\widehat{\theta }}_2=0.2045$	$\mathrm{BIA}=0.9462$
		${\widehat{\theta }}_3=-0.1916$	$\mathrm{MSM}=0.3714$
	Epanechnikov	$\widehat{h}=0.10,\widehat{p}=3$	$\mathrm{MFE}=0.0148$
		${\widehat{\theta }}_1=0.4011$	$\mathrm{MASE}=20.0997$
		${\widehat{\theta }}_2=0.2095$	$\mathrm{BIA}=0.9483$
		${\widehat{\theta }}_3=-0.1894$	$\mathrm{MSM}=0.3711$
	triweight	$\widehat{h}=0.11,\widehat{p}=3$	$\mathrm{MFE}=0.0157$
		${\widehat{\theta }}_1=0.3643$	$\mathrm{MASE}=19.7872$
		${\widehat{\theta }}_2=0.1853$	$\mathrm{BIA}=0.9457$
		${\widehat{\theta }}_3=-0.2175$	$\mathrm{MSM}=0.3954$
FATSM	Gaussian	$\widehat{p}=3$	$\mathrm{MFE}=0.011$
		${\widehat{h}}_1=1.5$	$\mathrm{MASE}=19.625$
		${\widehat{h}}_2=1.75$	$\mathrm{BIA}=0.964$
		${\widehat{h}}_3=0.05$	$\mathrm{MSM}=0.310$
	Epanechnikov	$\widehat{p}=3$	$\mathrm{MFE}=0.010$
		${\widehat{h}}_1=1.25$	$\mathrm{MASE}=17.625$
		${\widehat{h}}_2=1.73$	$\mathrm{BIA}=0.966$
		${\widehat{h}}_3=0.1$	$\mathrm{MSM}=0.33$
	triweight	$\widehat{p}=3$	$\mathrm{MFE}=0.0098$
		${\widehat{h}}_1=1.5$	$\mathrm{MASE}=16.745$
		${\widehat{h}}_2=1.77$	$\mathrm{BIA}=0.960$
		${\widehat{h}}_3=0.13$	$\mathrm{MSM}=0.345$

. The results reveal that the FNPTSM outperforms the FSPTSM and FATSM for the global land–ocean temperature dataset. Note that the best results of the proposed FNPTSM are obtained based on the Gaussian kernel and the best results of both the FSPTSM and FATSM are given when implementing the triweight kernel. The fuzzy data, along with the corresponding estimations related to the FNPTSM (based on the Gaussian kernel) as well as of the FSPTSM and FATSM (based on the triweight kernel), are visualized in Figure 3. In comparison to the FSPTSM and FATSM, the values predicted by the FNPTSM are closer to the fuzzy observations, which reveals that the proposed FNPTSM performs better for the global land–ocean temperature dataset.

5. Conclusions

Nonparametric statistical inference deals with situations where the functional relationships of the involved distribution functions are unspecified. In this regard, nonparametric time series models were broadly utilized to identify the “best fit” curve for a given time series of data. However, there are numerous situations where the available data are fuzzy rather than exact. In this paper, a nonparametric kernel-based time series model that relies on fuzzy data was introduced. The Nadaraya–Watson estimator was utilized to provide a fuzzy time series model within a three-stage procedure. Some popular goodness-of-fit measures have been implemented to investigate the performance of the fuzzy nonparametric time series model based on different kernel functions. The effectiveness and feasibility of the proposed time series model were also compared with the performance of existing time series models based on fuzzy data. Considering three common kernel functions (Gaussian, Epanechnikov, and triweight), the results indicated the superior performance of our proposed method in comparison to previous approaches. In addition to the performance aspect, the handling of the new nonparametric kernel-based time series model is much simpler than that of the previous methods, as we implemented an estimation procedure that is divided into three independent stages. In addition, our proposed time series model can be employed for arbitrary shapes of $LR$ fuzzy numbers. However, the model can be applied only for $LR$ fuzzy numbers, and thus it could be a promising future direction to develop a more general methodology that can handle arbitrary fuzzy numbers. Future studies could also focus on extending our approach to cases where the underlying time series data contain outliers. Finally, extending the proposed methodology for other nonlinear models such as wavelet-based or neural network-based time series models are further ideas for future research.