Next Article in Journal
Green STEM to Improve Mathematics Proficiency: ESA Mission Space Lab
Next Article in Special Issue
Relationships between Copper Futures Markets from the Perspective of Jump Diffusion
Previous Article in Journal
Machine Learning Approach for Modeling and Control of a Commercial Heliocentris FC50 PEM Fuel Cell System
Previous Article in Special Issue
Financial Development, Saving Rates, and International Economic Volatility: A Simple Model
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Chaos Analysis of the Dry Bulk Shipping Market

by
Lucía Inglada-Pérez
1,* and
Pablo Coto-Millán
2
1
Department of Statistics and Operational Research, Medicine Faculty, Complutense University, Plaza Ramón y Cajal, s/n Ciudad Universitaria, 28040 Madrid, Spain
2
Department of Economy, University of Cantabria, 39005 Santander, Spain
*
Author to whom correspondence should be addressed.
Mathematics 2021, 9(17), 2065; https://doi.org/10.3390/math9172065
Submission received: 23 July 2021 / Revised: 14 August 2021 / Accepted: 18 August 2021 / Published: 26 August 2021

Abstract

:
Finding low-dimensional chaos is a relevant issue as it could allow short-term reliable forecasting. However, the existence of chaos in shipping freight rates remains an open and outstanding matter as previous research used methodology that can produce misleading results. Using daily data, this paper aims to unveil the nonlinear dynamics of the Baltic Dry Index that has been proposed as a measure of the shipping rates for certain raw materials. We tested for the existence of nonlinearity and low-dimensional chaos. We have also examined the chaotic dynamics throughout three sub-sampling periods, which have been determined by the volatility pattern of the series. For this purpose, from a comprehensive view we apply several metric and topological techniques, including the most suitable methods for noisy time series analysis. The proposed methodology considers the characteristics of chaotic time series, such as nonlinearity, determinism, sensitivity to initial conditions, fractal dimension and recurrence. Although there is strong evidence of a nonlinear structure, a chaotic and, therefore, deterministic behavior cannot be assumed during the whole or the three periods considered. Our findings indicate that the generalized autoregressive conditional heteroscedastic (GARCH) model and exponential GARCH (EGARCH) model explain a significant part of the nonlinear structure that is found in the dry bulk shipping freight market.

1. Introduction

Since the 1960s chaos theory has played an increasingly important role in the development of many scientific fields such as economics and finance [1,2,3,4,5]. Interest in nonlinear modeling, and particularly in studying chaotic system, defined as a nonlinear deterministic process that appears to be random and presents sensitivity to initial conditions, has risen significantly. In fact, literature displays numerous publications searching for nonlinear and chaotic behavior in macroeconomic and financial time series [2]. Our research is included within this framework which aims to unveil the nonlinear dynamics in the shipping market. Specifically, we examine the existence of nonlinearity and chaos in dry bulk shipping rates time series.
Maritime transport plays an outstanding role in the worldwide economy, since over 90% of all world trade is transported by sea and it is the most profitable method of carrying goods and raw materials in masse throughout the world [6]. Specifically, ocean shipping freight rates constitute a significant proportion of the production inputs cost and determine to some extent the final commodities’ prices in overseas consumption centers [7]. Dry bulk freight rates are commonly represented by the Baltic Dry Index (BDI) [8] that has been proposed to be a measure of the shipping rates for certain raw materials [9]. The BDI is considered a global trade indicator [10] and serves as an indicator for the commodities, currency, and equity markets [8]. Moreover, it is considered an accurate “barometer” of economic activity and an efficient indicator for forecasting industrial production and economic growth because it captures changes in the early stages of the global supply chain. In summary, as Guan et al. [11] highlighted, BDI not only reflects the evolution of the dry bulk shipping market but can also mirror the state of the world economy and the trend in international trade.
Given their interest, it is not surprising that, as mentioned by Chen et al. [12], knowledge of freight rate dynamics has attracted a great deal of attention among researchers in the shipping community. Stopford [13] emphasizes the great importance of forecasting shipping prices for the stakeholders involved in the maritime industry: shipping companies, cargo owners, shipbuilders, bankers, or regulators. In addition, he points out that forecasts can achieve a high degree of complexity and in many cases, they have had limited success. Forecasting freight rates is a complex task because they are influenced by many features [14], such as political, economic, and natural conditions [11]. Among the many influences on the shipping market, Stopford [13] focused on factors affecting the demand for sea transport (i.e., world economy, seaborne commodity trade, and random shocks) and affecting the supply (i.e., world fleet, shipbuilding deliveries, and scrapping and freight revenues). Zeng et al. [15] suggested that the complexity as well as the non-stationary and nonlinear nature of freight rates series in the bulk shipping market escalate the difficulty in forecasting freight rates. In this sense, Chen et al. [12] mentioned that the dry bulk shipping industry is characterized as highly volatile. Moreover, Geman and Smith [16] focused on extreme volatility in the shipping markets, particularly in freight rates, because the inelasticity of supply to demand shocks when no inventory is available causes large fluctuations in the change through time. For the aforementioned reasons, accurate freight rate prediction presents an important challenge for researchers [15,17].
However only a limited number of studies that unveil the existence of nonlinearity and chaos can be found (e.g., [17,18,19]). Specifically, Goulielmos and Psifia [17] found a nonlinear behavior and low-dimensional chaos in the time charter freight rates series dynamics using the Hurst exponent, BDS test and Lyapunov exponent. Goulielmos and Psifia [18] detected a nonlinear behavior in the freight rates time series dynamic by applying the BDS test. Goulielmos et al. [19] found evidence of chaos, applying techniques like rescaled range analysis, correlation dimension and Lyapunov exponent. In conclusion, previous research on chaos has used methods from other sciences that are not appropriate for diagnosing chaos (for example BDS test or Hurst exponent) or that have not been adapted to take into account the specific characteristics of financial time series like noise and the limited and low sample size (e.g., correlation dimension method or Lyapunov exponent test). Hence, their conclusions could be misleading [20,21,22]. It can be concluded that the existence of chaos in ocean freight rates continues to be an open and important issue. Its relevance resides in that the existence of low-dimensional chaos could enable a reliable short-run prediction, but not in the long-run, since a chaotic system presents sensitivity on initial conditions and therefore is unstable [23].
From an economic point of view, the existence of chaos has several important implications. As Brooks [24] points out, these consequences affect aspects such as model appropriateness, market efficiency and predictability. LeBaron [25] argues that chaos in economics had a broad array of potential applications such as forecasting or understanding international business cycles. The presence of chaos can also provide a potential explanation for seemingly random fluctuations in the economy and financial markets [26] and a wider range of time series behavior. Specifically, linear models may not adequately capture the sharp movements and large fluctuations that characterize financial and economic series and hence are not appropriate to describe economic phenomena like stock market bubbles and crashes, depressions, and the occurrence of business cycles. On the contrary, chaotic models may be well suited to explain these behavioral patterns [27]. Another important implication of the presence of chaos is that controllability of the system is feasible, since there is some deterministic structure underlying the data [28]. Interestingly, chaos also represents a radical change of view on business cycles [28]. The traditional exogenous approach to economic fluctuations assumes that economic equilibria are determinate and intrinsically stable, so that in the absence of persistent exogenous shocks the economy will tend to a steady state, but due to stochastic shocks a stationary pattern of fluctuations is observed. In contrast, with the new approach based on chaos theory, business cycles are endogenously explained and are based on the marked nonlinear deterministic structure that can characterize the economic system.
Concerning the degree of predictability of the series, if the system is completely random, its behavior is not predictable at all, while if it exhibits chaotic behavior, we can forecast the system in short periods of time [29]. However, when chaos exists, although short-run forecasting may be more accurate, long-run forecasting then becomes impossible due to the sensitivity to initial conditions. In addition, traditional linear forecasting methods would not be successful and nonlinear models should be applied [26]. The presence of chaos also indicates the possibility of using nonlinear trading rules and is likewise understood as a hallmark of an inefficient market [30]. In practical terms, this implies that in financial markets accurate short-run forecasts can be achieved and investors can make profits with an adequate strategy, while the long-run behavior of these time series is not predictable in any way [31].
As noted above, linear models have proven to be ill-suited for describing and fitting many economic phenomena such as the business cycles. In contrast, the nonlinear approach can capture the features of economic and financial series and their sudden fluctuations, and hence plays a relevant role in economic modeling [32]. Likewise, dynamical systems involving time delays have applications in a number of fields such as biology, population dynamics and economics [33]. The delayed model can produce oscillating trajectories for the solutions under certain conditions and explain the business cycles [34,35]. Delay differential equations have derivatives depending explicitly on the value of variables at times in the past and are adopted to represent systems with time delay. In general, delay differential equation models can be more accurate than those based on ordinary differential equations when it is necessary to capture oscillatory dynamics. In the last few years, the literature on the oscillation theory of delay differential equations is growing very quickly and among the topics studied the oscillation of solutions has been the most extensively investigated [36]. The core of the theory lies in setting the conditions for the existence or non-existence of oscillatory solutions. This is an area of research with interesting applications in many fields and closely related to chaos, which is characterized by irregular and unpredictable oscillations [37]. Hence, oscillatory behavior of the solutions of a class of nonlinear delay differential equations is an important area for future research.
The basic aim of this research is to investigate the underlying dynamics of ocean freight rates as measured by the BDI. By using daily closing values from January 1990 to July 2013, we tested for the existence of nonlinearity and low-dimensional chaos. The paper also examines other issues such as the existence and modelling of volatility clustering. The proposed methodology considers the characteristics of chaotic time series, such as nonlinearity, determinism, sensitivity to initial conditions, fractal dimension and recurrence. Each characteristic is tested by the corresponding method. The existence of nonlinearities is tested by the application of a broad set of methods (Keenan, Tsay, White and BDS). In the same way, four methods (correlation dimension, Matilla-Garcia and Ruiz Marin (MGRM) test [38], Lyapunov exponent test, adapted for noisy series [21], and recurrence plots) have been applied to test for chaos. To make the conclusions of this work more robust, besides analyzing the series over the whole period, we have also examined the chaotic behavior over three sub-periods, which were selected by the evolution of volatility.
To this end the following methodology is applied to all time series [32]. First, to evaluate stationarity of the BDI series, unit root tests were performed. Once it had been established that BDI series was I(1), the analysis was carried out on daily returns computed as first difference of the logarithmic price level. Linear dependence was then filtered out of the data by applying autoregressive moving average (ARMA) filters based on the Box–Jenkins methodology [39]. Next, we applied a broad array of procedures to the residuals obtained to detect the existence of nonlinearity. If nonlinear dependence was detected, since linear structures had already been removed using the best fit ARMA model, it was indicative of some type of nonlinear dependence in the returns series resulting from a nonlinear stochastic system or a nonlinear deterministic chaotic system. Stochastic nonlinearity could be due to the presence of a volatility cluster, i.e., large variations tend to group together, which is a common characteristic of financial time series. In such a case, both generalized autoregressive conditional heteroskedasticity (GARCH) [40] and exponential (EGARCH) [41] models were fitted. Afterwards, as a distinctive feature of proposed methodology against published methods, the existence of chaotic structure was examined by a series of robust procedures, including those that are ideally suitable for noisy time series analysis, e.g., Lyapunov exponent test adapted for noisy series, MGRM test based on permutation entropy, and recurrence plots.
Furthermore, the novelty and main contribution of this study lies in the fact that, as far as we know, it is the first time that these robust methods have been applied to detect chaos to such a relevant variable in the world economy, the freight shipping rate. The findings suggest that nonlinearity is present in the underlying process of the BDI although unlike previous studies, we found no sign of chaotic structure. In terms of the benefits of our research for both theory and practice, our results make a significant contribution to shedding light on the underlying behavior and dynamics in shipping rates. In general, our insights should be of great interest to all stakeholders concerned with the shipping industry, including, among others, stevedores, charterers, shipowners, brokers and investment banks.
The rest of the paper is organized as follows. The next section describes the methodological framework applied. The penultimate section presents the data sources used and the main results. Finally, in the final section, we present some conclusions and suggest future lines of research.

2. Methodology

2.1. Methodological Framework

The methodological framework from other research [32,42] has been adopted to study the presence of nonlinearity and chaos. The methods employed herein are presented in Table 1, and their main characteristics are described below.
Firstly, we define linear and nonlinear series as elsewhere [32]. A stochastic process ( X t ,   t T ) is said to be a linear process if for every t T   X t = j = 0 β j ϵ t j where a 0 = 1, ϵ t ,   t T is a process of independent and identically-distributed random variables (iid) with E X t = 0 ,   E X t 2 = σ 2 and j = 0 β j < . If a process does not satisfy this condition, then is said to be nonlinear.
The methodological framework includes the following steps.
(i)
Stationarity testing
Stationarity was studied by means of three different techniques (Table 1). Specifically, the Augmented Dickey Fuller test (ADF) [43], Phillips–Perron test (PP) [44] and Kwiatkowski–Phillips–Schmidt–Shin test (KPSS) [45] are applied to time series.
In first instance, once it was tested that each of the four series considered was not stationary, we obtained the log-returns (from now on returns) from each series defined as ln R t = ln X t ln X t 1 such that X t is the time series data at time t. This new series called “return” has important properties such as stationarity. This property is essential, since it is a prerequisite in a large percentage of modelling techniques and in the methods that will be used for the analysis of nonlinearity or chaotic behavior.
(ii)
Linear modelling
Next, ARMA filters were applied to remove the linear dependence present in the return series. The specific model is chosen as the ARMA (p,q) model which presents the lowest value according to the Schwarz information criterion [53]. Following the process described by Barkoulas et al. [54], the residuals were studied. If the residuals were correlated, the order of the ARMA model was increased, by increasing the magnitude of p and q. Once the optimal linear model was determined, the residuals were extracted from it. The new residuals were called the ARMA series.
(iii)
Nonlinearity testing
We used several methods to supplement the limitations of each one. We applied (Table 1): Keenan’s test [46], Tsay’s test [47], White’s test [48] and BDS [49,50].
Before describing the main characteristics of the methods employed in this research, it is worth explaining the process of reconstructing the state space, as well as the concept of attractor.
Following [54], the concept of attractor is defined as follows. Given a discrete dynamic system of the type:
x t = F x t 1 , x R n
such that P an open subset of R n and F : P R n is a function [55]; a closed invariant set A P is an attracting limit set of P if there is an open neighborhood Q of A , such that the limit set of iterates is A , x ϵ Q when t .
Empirically, a time series of scalar observations x t , representing the multidimensional system of Equation (1), is usually observed. Using Takens’ embedding theorem [56] it is feasible to recover the dynamics of the system (i.e., the original trajectory) from the observed x t time series. An m-dimensional vector is defined from x t , and constructed as follows:
X t m = x t , . , x t + m 1 = g x t . , g F m 1 x t J m
where F m 1 means that F is applied a total number of m − 1 times to itself. Thus, our goal is to reconstruct the state dimension space by expanding the one-dimensional signal x t into an m-dimensional phase space such that each observation in the signal x t is substituted by the corresponding vector X t m in Equation (2). Takens’ embedding theorem states that for each pair F , g the map J m :   R n R m will be an embedding for m 2 n + 1 . If the embedding dimension m is of sufficiently large size with respect to the attractor dimension, Takens’ theorem ensures that a diffeomorphism exists between the original and the reconstructed attractor. In summary, according to the above theorem it can be concluded that both attractors represent the same dynamical system in different coordinate systems if m 2 n + 1 .
Nonlinearity tests used in this research are well known in the literature. Thus, only the main features of the BDS method, which has statistical power for detecting nonlinearity and other types of data dependence [21], are presented here. Interestingly, the BDS test can be used to test for nonlinearity if linear structure is removed from the series by fitting an ARMA model. Indeed, if the test is performed on the residuals from a linear model (ARMA process), rejection of the null hypothesis (iid) will denote the presence of linear dependence in the series. The BDS test is constructed as follows [49,50]. Consider a time series x t t = 1 T . According to Takens’ theorem it is possible to reconstruct the state space from x t . The method of delays is used for reconstruction. Thus, the reconstructed state vector or m-history is   X t m = x t , x t + τ , . , x t + m 1 τ where τ is the lag or reconstruction delay, and m is the embedding dimension of the reconstructed space. The reconstructed trajectory, X, of N m-dimensional vectors X t m can be expressed as a N x m matrix with N = T m 1 τ .
Under the null hypothesis, the BDS statistic is derived using the correlation integral, which measures how many times the temporal patterns in the data are recurrent.
V T , m , ε = N C T , m , ε C T , 1 , ε m σ ^ T , m , ε
being C T , m , ε the correlation integral which can be stated as:
C T , m , ε = 2 N N 1 t < s H ϵ | | X t m X s m | |  
where ‖.‖, m and H represent a norm operator, the embedding dimension and the Heaviside function, respectively.
The correlation integral asymptotically follows a standard normal distribution. σ ^ T , m , ε is the standard sample deviation of C T , m , ε C T , 1 , ε m . Assuming a time series as iid, BDS tests for the null hypothesis that C T , m , ε = C T , 1 , ε m , which is tantamount to the null hypothesis of iid versus an undefined alternative [1].
(iv)
Volatility modelling
In the case of identifying the existence of nonlinear dependence in the residual series from ARMA models, it could be due to the existence of a volatility cluster. To corroborate this hypothesis, conditional variance is modeled by fitting ARCH-type models (GARCH and EGARCH). The selection of the best model is carried out using the same method that was employed for the ARMA model. Once both types of model have been estimated, their standardized residuals are obtained from their conditional standard deviations (hereafter GARCH and EGARCH series). The techniques for detecting the existence of nonlinearity mentioned above are applied to these new GARCH and EGARCH series to contrast whether all the nonlinear dependence has been removed by the process performed.

2.2. Tools to Detect Chaotic Regime

Below, we describe the main characteristics of the tools used for detecting a chaotic regime.

2.2.1. Correlation Dimension

Fractal dimension represents a measure for the level of complexity in the chaotic attractor [1]. This measure has the property of remaining unaltered by the state space reconstruction process. A characteristic of chaotic systems lies in the fact that the fractal dimension is not an integer value [1]. Therefore, estimation of this parameter can be used as a test to detect the presence of chaos. Among the various algorithms for estimating the fractal dimension, the Grassberger and Procaccia algorithm [51] for obtaining the correlation dimension has practical advantages due to its ease of implementation. In fact, with this correlation dimension method, we can distinguish stochastic series from chaotic series [29].
Analogously to the BDS test, this method is also based on the reconstruction of the state space from the observed series and on the correlation integral. Briefly, let x t t = 1 T be a time series and the sequence of N = T m + 1 m-dimensional vectors or m-histories derived from it, X t m = x t , , x t + m 1 that yields the reconstructed series. The Grassberger and Procaccia algorithm [51] is employed to compute the dimension of the reconstructed attractor. The method relies on the correlation integral which is defined as follows:
C ϵ = 2 N N 1 t < s H ϵ | | X t m X s m | |  
where m, ‖.‖ and H represent, respectively, the embedding dimension, norm operator and the Heaviside function.
To determine the correlation dimension, it is necessary to determine how C   ε varies as ϵ varies. As the value ϵ rises, the value of C   ε increases due to the greater number of nearby points to include. Grassberger and Procaccia [51] proved that in the case where ε is small enough, C   ε may be well-approximated by C   ε ~ ε v . In other words, when ε→0, C   ε grows at rate v . The estimate of v when m , provides the correlation dimension (CD). In practice, the CD of a dynamic system is determined by estimating the slope of the regression of ln   C   ε versus log ϵ and an intercept for small values of ε , and depends on the chosen embedding dimension.
In the case where the data are purely stochastic, the value of the correlation dimension will be equal to m for all values of m. If the series is deterministic in nature, the estimated slope will plateau at a point, not rising as m rises. This level of “saturation” of the slope is precisely equal to the estimated correlation dimension for the underlying process.

2.2.2. Lyapunov Exponent

Sensitivity to initial conditions is widely recognized as a core feature of chaos. A checkable quantitative criterion for the sensitivity to initial conditions can be given by the so-called Lyapunov exponents [57]. They calculate the rate of divergence of closely spaced trajectories whose initial conditions only diverge in an infinitesimal amount. LE is one of the most widely used methods to investigate the existence of chaos in time series [23]. The largest Lyapunov exponent can be employed to gauge the separating rate of nearby trajectories and estimate the level of chaotic structure that is present in a nonlinear dynamic system [58]. Under certain conditions, the sign of the coefficient or exponent is a signal of the presence of chaos in the studied series [57]. Specifically, our attention is restricted to systems with bounded metrics, or closed, bounded invariant subsets of unbounded systems. For systems with bounded fluctuations, if the Lyapunov exponent is positive, the orbit exhibits sensitive dependence on initial conditions, and this is a necessary condition of chaotic behaviour [59]. In summary, as Nychka et al. [60] and Bask et al. [61] pointed out, a bounded system with a positive Lyapunov exponent is an operational definition of chaotic behavior in nonlinear dynamic systems that appear to have random and unpredictable behavior.
There are many algorithms in the literature for the estimation of this coefficient. The algorithm proposed by BenSaïda and Litimi [21] is followed in this research. Specifically, it has been empirically proven that this methodology to determine the Lyapunov exponent (λ) is the one that behaves best for noisy series [62]. Indeed, literature shows that there are two main types of technique, which are founded on the reconstruction of the spatial state by means of delay coordinate embedding approaches, to estimate λ from experimental data. Direct methods are based on the computation of the divergence growth rate between two trajectories with an infinitesimal difference in their initial conditions [63]. One of the relevant limitations of this method is that it is sensitive to both sample size and level of noise in the data [64]. By contrast, the Jacobian-based approach can provide robust estimates of Lyapunov exponents despite the existence of noise [65,66]. The Jacobian method relies on non-parametric regression to estimating the Jacobians and λ. However, as for a scalar time series, the map that triggers the process is generally not known; the Jacobian matrix cannot be obtained and, therefore, the Lyapunov exponent cannot be determined. To do so, it is necessary to approach the unknown chaotic map by means of a given function. McCaffrey et al. [65] compared the results of various approaches: radial basis functions, neural networks, thin-plate splines, and projection search. On the basis of their simulations, neural network behaved as the best regression method for the case of chaotic systems with noise. In this line, Bailey et al. [66] suggested a regression method that involves using neural networks. Simulation results for a Henon system containing noise showed that neural network regression method provides accurate estimates for the Lyapunov exponent.
Following [21], the methodology comprises the following steps: consider a time series x t t = 1 T . A noisy chaotic system can be expressed as follows:
X t = f x t τ , x t 2 τ , ,   x t m τ + ε t
where t is the time script and m, τ, ε t and f represent the embedding dimension, the time delay, noise added to the series and an unknown chaotic map, respectively. The Lyapunov exponent (LE) is set as:
λ ^ = 1 2 M ln v 1
where M T is the “block-length” or number of equally spaced evaluation points that are used to calculate the Lyapunov exponent and v 1 is the largest eigenvalue of the matrix T M U 0 T M U 0 , being T M = t = 1 M 1 J M t ,   J t the Jacobian matrix of the chaotic map and U 0 = 1 ,   0 ,   ,   0 . Since we usually do not know the chaotic map F, the Jacobian matrix could not be estimated. Therefore, it is necessary to approximate the unknown chaotic map by the appropriate function.
As in [18], a neural network function, with a q-dimensional hidden layer and a hyperbolic tangent function as the activation function, is adopted. It is estimated by the nonlinear least squares method, varying m from 1 to 8, in order to subsequently compute the LE spectrum. Thus, the chaotic process is estimated by using the equation below:
x t α 0 + j = 1 q α j tan h β o , j + i = 1 m β i , j x t i τ + ε t
The triplet (τ, m, q) defines the complexity associated with the process and is selected as the one that determines the largest value for λ. Then, from the asymptotic distribution of λ, the test for the presence of chaos is constructed [63].

2.2.3. Matilla-Garcia and Ruiz Marin (MGRM) Test

The MGRM Test [38] is a test for the deterministic process as opposed to the stochastic one which is rooted in symbolic dynamics and entropy. In particular, the concept of permutation entropy, which has its origin in symbolic dynamics, is the foundation of the MGRM test. The basic idea of symbolic dynamics consists of partitioning the phase space into a finite set of regions and to assign an alphabetic symbol to each region. A relevant property of symbolic dynamics is that essential features of the underlying dynamics, namely deterministic or stochastic nature or complexity, are preserved. Likewise, entropy provides an insight into the unpredictability of the system under study, which is a key feature of complex systems.
In the first instance the notation is introduced: let x t t T be a stationary time series, x t t I an observation with I = 1 , , T and m the embedding dimension with m 2 . We will define the ordinal patterns for “m”. As we already know, the scalar time series can be embedded in an m-dimensional space as follows: x m t = x t , x t + 1 , , x t + m 1   for   t T . For a given time value t, the ordinal pattern corresponding to embedding dimension m is set to be the unique permutation π m t r 0   r 1 r m 1 of the set 0 ,   1 ,   .   .   . ,   m 1 that satisfies:
x t + r 0 x t + r 1 x t + r m 1  
r s 1 < r s   if   x = x t + r s  
By means of Equation (9) it is guaranteed that the permutation defined by Equation (10) is unique. Therefore, the vector or m-history X m t is transformed into a unique symbol π m t . In fact, π m t shows us how the ordering of the times: t + 0 < t + 1 < < t + m 1 becomes the order of the correspondent values being analyzed. The core point is to split the state space in which the dynamical process unfolds into a finite number of sets exploiting the time dependent information held in the m-history X m t R m . According to the above definition, the partitions depend on how the ordinal structure corresponds to the m-history is. In particular, π m t = π m s ,   s t , if and only if for all k , I   ϵ   0 , 1 , , m 1 with k I it holds that x t + 1 x t + k x s + 1 x s + k .
In general, given a time series x t t T , all m! permutations of order m are taken as potential order types of m distinct numbers. Then, for an embedding dimension parameter m, the relative frequency or unconditional success probability p π of each symbol or permutation π can be set as:
p π = c a r d t   |   0     t     T m 1 ,   Y m t   h a s   t y p e   π   T m + 1
Given an embedding dimension   m 2 , Shannon entropy that stands for the m! distinct symbols, is defined as follows:
h m = i = 1 m ! p π i log p π i  
The test is constructed as follows: Let m be embedding dimension and T be the number of observations and fix w , k such that w = m ! k . Next, the subsets W j that must verify W 1 W 2 W n in the following way W j = W j 1   w   symbols   chosen   at   random   in   S m   \   W j 1 for j = 2, …, k are constructed.
Subsequently the next permutation entropy’s function is calculated as:
h W j m = π j w j p π log p π
When analyzing the permutation entropy’s function values, it is possible to distinguish and identify deterministic systems. In this process no additional information is obtained by an increase in the number of symbols considered. In contrast, in non-deterministic process, this information or complexity is increased.
This latter property is tested in the following way: Let d h W j m be the h W j m slope:
d h W j m = h W j + 1 m h W j m log j + 1 j
When considering random process, the numerical slope of permutation entropy will rise with log(jw), while this will not be the case for chaotic or regular processes.
The property described above that identifies deterministic systems h W j m is tested by carrying out the following regression:
d h W j m = α 0 + α 1 j + ε j   f o r   j = 1 , 2 , ,   k 1
where ε j is white noise. As a result, the estimated parameter α ^ 1 can be employed to assess how d h W j m scales with j. In mathematic notation, the test is the following one:
H 0 α 1 = 0 H 1 α 1 > 0
Indeed, regression (15) can be considered as a simple symbol-trend model. The same as for the basic time-trend model, the ordinary least squares estimate α ^ 1 is such that the standard t-test of H0: α ^ 1 = 0 is asymptotically valid. If the coefficient obtained is null, then the series is deterministic; otherwise (greater than 0), it is stochastic [38].

2.2.4. Recurrence Plots

A recurrence plot (RP) reveals the existence of recurrent and intermittent patterns in time series. First proposed by Ekmann et al. [52,67], it has been widely applied in the characterization of dynamic systems. This topological method shows the hidden structures of time series from a qualitative perspective. The plots are constructed by assuming mutual distances that belong to the same path in the reconstructed state space [42]. The construction is done as follows.
Let m be the embedding dimension and X t m the vector m-dimensional in the reconstructed state space at time t ranging from 1 to k. The recurrence matrix is generated by comparing each embedded vector X i m with the other X j m . A point is drawn if this comparison is less than a value ε for a specific distance. That is, if the condition is met: | | X i m X j m | | < ε . Recurrence of a state between t = i and t = j occurs when X j m is close enough to X i m . The RP can, therefore, show which vectors are close together and which are far apart.
The diagonal structures of the RP define the extent to which one part of the path is relatively near to another at a distinct time. For the case of a deterministic chaos system, short lines parallel to the main diagonal are visible. However, in stochastic systems, short line segments are missing, and uniformly distributed points are displayed.
This technique is independent of some constraints such as sample size, noise, and stationarity [1]. It also provides additional information about the structure of the attractor, since the plot preserves the temporal order of the series, allowing the place and periodicity of the periodic orbits to be known.

3. Empirical Results

Applying the methodological framework described above, the main results obtained are shown below.

3.1. Data

The data set corresponds to daily observations of the BDI comprising the period between 1 January 1990 until 1 July 2013 and was obtained from the Bloomberg (http://www.bloomberg.com/ (accessed on 26 May 2014) and Baltic Exchange Web sites (http://www.balticexchange.com/ (accessed on 26 May 2014)). In addition, as a complement to the above and after checking that the initial series had periods of different volatility, it was decided to divide the original series into three new series, taking into account the following sub-periods: (i) BDI1 series (from 1 January 1990 until 31 December 2002); (ii) BDI2 series (from 1 January 2003 until 31 December 2007) and (iii) BDI3 series (from 1 January 2008 until 1 July 2013).
Figure 1 shows that the BDI remained relatively stable, following a course lacking in significant shocks until 2003. Then, there was a dramatic increase in volatility, adopting a behavior that mimics financial bubbles. After the ocean freight boom (2005–2006 period), on 20 May 2008, the BDI attained its maximum value (11,793) since its introduction in 1985. However, it immediately decreased 94% during the period of May to December 2008, in the midst of the complete turmoil of the financial crisis. It was not until 5 December 2008, that it registered a value of 663 points, corresponding to its minimum value since 1985. This decrease is partially related to the reduction in demand for global goods that translated into a fall in the demand for transportation and coincides with an increase in the supply due to the deliveries of new ships ordered during the boom period of the global economy.
Table 2 presents the main descriptive statistics of the four series studied. All series show a positive skewness coefficient. Moreover, it can be concluded that the data do not come from a normal distribution, as reflected by the Jarque–Bera test (p-value < 0.05).

3.2. Returns

To evaluate the stationarity of the BDI series, unit root tests are performed [68,69]. Specifically, we examine the stationary property by applying the Augmented Dickey–Fuller test [43], the Phillips–Perron test [44] and the Kwiatkowski–Phillips–Schmidt–Shin contrast [45]. After testing that the four series are not stationary (see Table 3), the returns from each series were obtained by considering the logarithms of the ratio of two consecutives prices. Figure 2 shows the time development of the returns. The peak in 1999 matches with the methodological change implemented in the index definition. Return series shows a negative skewness coefficient, that is, the distribution is left-skewed, implying that upward jumps of larger size are less frequent than downward jumps of smaller size, as do financial time series [70]. It is worth noting the high value of the kurtosis coefficient, which implies that large jumps are more likely to happen more frequently than in the normal distribution. This result is not altogether surprising because the BDI series presents similar characteristics to financial series, as previously stated [71].
In the evolution of the return series corresponding to each of the three sub-periods, two peaks corresponding to the 1997 Asian crisis (BDI1 returns series) and the 2008 world recession (BDI3 returns series) are displayed. The three series present a positive skewness coefficient and an excess of kurtosis, with the latter reaching a greater magnitude in the third series (BDI3 returns series). According to the results from the Augmented Dickey–Fuller, Philips–Perron and KPSS tests, all the three returns series are stationary (Table 3).

3.3. Linear Dependence

As some methodology (for example the BDS test) is not robust for the presence of linear relationships, the linear dependence was removed using ARMA models. Following the criterion described by [54], the best ARMA model was selected for each series (Table 4). The selected fitted models for each of the series (ARMA(9,0); ARMA(5,0); ARMA(3,0) and ARMA(7,0)), for the series BDI, BDI1, BDI2 and BDI3 respectively, are shown in Table 4 that displays all the models fitted in this research. Details of each model are presented in Appendix A (Tables A1, A4, A7 and A10). According to the results from the Augmented Dickey–Fuller, Philips–Perron and KPSS tests applied on the residuals of the selected model, and shown in Table 3, all the residuals series (ARMA) are stationary.

3.4. Nonlinear Dependence

Next, Keenan, Tsay, White and BDS methods were used to study nonlinearity. The principal advantage of using a broad set of instruments is to obtain the most information possible about the nature of the series, since to date none of the methods has proved successful in detecting all types of dependence. Nevertheless, some tests, such as the BDS or White are more efficient in detecting dependency [3]. Results by applying Keenan, Tsay and White methods are exhibited in Appendix B (Table A13 and Table A14).
The BDS test was applied once the linear dependence was eliminated, as previously suggested by [73]. In this way the BDS test served as an indirect method of analyzing nonlinearity: if it rejected the null hypothesis, then the series was proven to be nonlinear. This procedure is theoretically feasible, since the calculations on the residuals of an autoregressive model do not lose the relevant information derived from the original series if the latter comes from a nonlinear chaotic system [73]. The BDS test is applied taking into account different embedding dimensions (range: 2–9). The results are presented in Appendix B (Table A15 and Table A26).
Table 5 presents a summary of results for the nonlinear analysis. Results from all procedures suggest the existence of nonlinearity in all the four returns series. Likewise, if the results of the tests with the highest power are considered (White and BDS), all ARMA series seem to present some type of dependence. In particular, the BDS test shows significant results in all ARMA series for the different embedding dimensions and ϵ values being considered.

3.5. Volatility Clustering

The next step consists of modelling the conditional variance by fitting ARCH family models (GARCH and EGARCH). In first instance, the best GARCH (p,q) model was selected so as to adhere to the previously described criterion [54]. After fitting each model, standardized residuals were obtained (from now on GARCH and EGARCH series). Unlike GARCH models, EGARCH models can be used to estimate conditional variance, taking into consideration the sign of the past period’s innovation. The asymmetric response in the conditional variance is successfully captured by this class of model and, therefore, they are well suited candidates for modeling financial series. The chosen models were EGARCH (2,3) and EGARCH (3,1) for BDI and BDI1 series, respectively, and EGARCH (1,1) for both BDI2 and BDI3 series (Table 4). Details of each model are presented in Appendix A (Table A2, Table A3, Table A4, Table A5, Table A6, Table A7, Table A8, Table A9, Table A11 and Table A12). The standardized residuals of all the series are less leptokurtic than those obtained from ARMA models (Table 2).
Since nonlinearity is a necessary, but not sufficient, condition for chaotic behavior, its existence was first analyzed in residuals of the volatility models (Table 5 and Table A13, Table A14, Table A15, Table A16, Table A17, Table A18, Table A19, Table A20, Table A21, Table A22, Table A23, Table A24, Table A25 and Table A26 in Appendix B). In general, linear hypothesis is not rejected, but nonlinearity was found in some cases. Indeed, both the White and BDS tests show some remaining dependence in the BDI series, which is compatible with the existence of chaos in the EGARCH model residuals. Almost all tests do not reject linearity for the BDI1 series in both the GARCH and EGARCH models. Likewise, almost all tests do not reject linearity for the BDI2 series in both the GARCH and EGARCH models. Finally, for the BDI3 series, both the White and BDS tests show some remaining dependence, which is compatible with the existence of chaos in both the GARCH and EGARCH model residuals.

3.6. Chaotic Behavior Analysis

The following methods were used to test the presence of chaotic behavior. Results are presented in Appendix C.

3.6.1. Correlation Dimension

The correlation dimension (CD) [51] quantifies the degree of complexity of a system and distinguishes a deterministic system from a stochastic one. The CD has been calculated for both the residuals from linear filters and the residuals from nonlinear models, considering different dimensions of embedding (m = 1, 2, …, 8). Likewise, for its analysis, the optimal value of the delay time has been considered, according to the criterion of minimum mutual information [74].
The results for the BDI series are shown in Appendix C (Table A27) and suggest that the correlation dimension generally grows as the embedding dimension rises. However, values of the CD, both for the return series and for the ARMA(9,0) and GARCH(1,1) series, do not converge to a fixed value, suggesting the absence of saturation, even though the increase of the CD is less than one. In the case of the EGARCH(2,3) series, this increase is greater than one, suggesting a stochastic behavior of the series. In the three subseries (BDI1, BDI2 and BDI3) analyzed, a similar result is observed: CD increases as the embedding dimension increases, being this increase greater for the case of nonlinear models. In short, the results suggest that in all cases there is sufficient evidence against the existence of a strange attractor.

3.6.2. Lyapunov Exponent

Lyapunov exponents (LE) are used to determine the sensitivity to initial conditions, a feature of chaotic processes. They calculate the rate of divergence of closely spaced trajectories whose initial conditions only diverge in an infinitesimal amount. Under certain conditions, a system is defined as chaotic when the associated LE is positive.
We calculated the maximum LE as described by [21] because this technique is able to deal with noisy data. Although we considered different values for the parameters of the neural network, negative significant Lyapunov exponents were obtained in all cases (Table A28 in Appendix C). Hence, we do not find evidence for the presence of chaotic structure. For all the series studied, the null hypothesis of existence of chaotic dynamic is strongly rejected at a significance level of 5% (p-values < 0.000).

3.6.3. MGRM Test

Subsequently, the MGRM test [38] has been applied to the series obtained in the nonlinear models. This test does not need to reconstruct the attractor of the system as a previous step and the analysis is done directly, using the data. The test is constructed as follows. First, the Shannon entropy values are calculated. Then, the information or complexity is estimated by a linear regression that used the entropies as the explanatory variable. The coefficient from the entropies is then analyzed, as in deterministic series the complexity derived from the entropy of permutation does not increase when the number of symbols increases, once the saturation point is reached. To conduct the test, the parameters were fixed such that m = 4, k = 2, and w = 12. The results were positive in all the cases (Table A29). Thus, the hypothesis of determinism is rejected.

3.6.4. Recurrence Plots

Visual recurrence analysis is based on Eckmann’s [52] definition of a recurrence plot. The degree of complexity and the existence of chaos is analyzed by generating the recurrence plots. This kind of analysis can capture the recurrence property of states, one of the essential properties observed in chaotic systems. The points above the main diagonal, representing the distance between the same vector and each embedded vector in the phase space and, therefore, segments parallel to the diagonal, would indicate a chaotic behavior [54].
Figure 3, Figure 4, Figure 5 and Figure 6 show the recurrence plots for the four series considered (BDI, BDI1, BDI2 and BDI3). For their elaboration, the optimal values of the two parameters, m—the embedding dimension and τ—the considered delay, as previously described, have been considered. They were estimated using the mutual information criterion and the false neighbours method, respectively [55,74].
The recurrence plots (RPs) for the different series considered are shown in Figure 3, Figure 4, Figure 5 and Figure 6, for the BDI, BDI1, BDI2 and BDI3 series, respectively. For its elaboration, the Euclidean distance has been chosen and the radius cut-off point (ε) is set as 10% of the largest distance across all points in the reconstructed state space. This value is adopted in several works, for example in [54].
The RPs corresponding to BDI series (Figure 3) do not exhibit any lines parallel to the main diagonal, which could suggest chaotic behavior [54]. Moreover, more intense horizontal and vertical lines are found in the uppermost right corner, indicating dynamic alteration periods, most likely due to the last worldwide economic crisis. In the nonlinear models a randomly distributed pattern of points that is characteristic of stochastic series is observed.
All RPs corresponding to BDI1 series (Figure 4) show a string of vertical and horizontal lines, located slightly to the right of the middle of the sample period, coinciding with the 1997 Asian crisis. Neither is there any indication of a chaotic component, since the previously described characteristic is not observed. The general pattern of all RPs corresponding to BDI2 series (Figure 5) is similar to the previous ones, with no lines parallel to the main diagonal. However, in this case, a greater number of horizontal and vertical lines are observed than in the previous cases. All RPs of BDI3 series (Figure 6) display the existence of time periods characterized by a prolonged dynamic alteration, especially at the beginning of the sample space and which is probably due to the last worldwide financial crisis.

4. Conclusions

This paper analyzes the underlying dynamics of dry bulk shipping rates using daily data from the Baltic Dry Index. The entire sample has been divided into three periods, chosen according to the volatility behavior of the series, to make the conclusions of this empirical work more robust. Specifically, we tested for the existence of a chaotic regime in the shipping market by applying metric, topological and permutation entropy-based approaches. From a comprehensive view a great number of methods have been used (four for nonlinear analysis and four for chaotic behavior analysis) including the most suitable for noisy time series. The adopted methodology is based on each method checking one of the properties that characterize chaotic behavior. Our results support the existence of nonlinearity, which is not consistent with chaos, as measured by the four indications of chaotic behavior, throughout all investigated periods. Specifically, the underlying system is of high dimensionality does not present sensitive dependence on initial conditions, as well as not being characterized by recurrent states and by deterministic nonlinear structures. In addition, we show that GARCH/EGARCH models explain a significant part of the nonlinear structure that is found in the dry bulk shipping freight market. Our findings are in accordance with the conclusions of other researchers who use novel and suitable procedures for noisy series to detect chaos in other financial time series (e.g., [21,54]). On the other hand, our conclusion contradicts the findings of previous studies that found evidence of low dimensional chaos in shipping market but using a number of techniques that are not suitable for detecting chaos (e.g., the Hurst coefficient) or that do not account for noise (e.g., [19]). Thus, their findings might be biased.
Regarding the benefits of our research for both theory and practice, results significantly contribute to illustrate the behavior and underlying dynamics of the maritime transport price. In general, our findings should be of great relevance to all stakeholders concerned with the maritime sector, such as institutional agencies, investors, charterers, shipowners, stevedores, and brokers. Important managerial insights can be drawn from analytical results. In particular, knowledge and forecasting of freight rates dynamics have a relevant role to improve the formulation of financial strategies, such as the pricing of derivatives and hedging instruments. Likewise, due to the volatile nature of freight rates to improve our volatility, forecasting accuracy is a major feature in successful risk management. Finally, since the BDI not only informs about the evolution of dry cargo shipping market but can also mirror the trends in international trade and consequently, the state of the world economy [11], it could be considered an accurate “barometer” of the economic activity and an efficient indicator for forecasting industrial production and economic growth. Hence, our findings are useful for the design of optimal macroeconomic policies.
Further research on the nature and forecasting of shipping rates is still an outstanding issue, as a limited number of studies are available. Future research is needed to examine whether other types of nonlinear model, both parametric and non-parametric, might be able to account for the remaining nonlinear dependence in the series. Also, studying the existence of nonlinearity and chaos in the volatility series would help to complete this work. Finally, a new investigation could consist of studying the effects of the COVID−19 outbreak on the dynamics of the analyzed series, as well as specifically investigating how this phenomenon could have modified the results of this research.

Author Contributions

Conceptualization, L.I.-P. and P.C.-M.; methodology, L.I.-P. and P.C.-M.; software, L.I.-P.; validation, L.I.-P. and P.C.-M.; formal analysis, L.I.-P. and P.C.-M.; investigation, L.I.-P. and P.C.-M.; resources, L.I.-P. and P.C.-M.; data curation, L.I.-P.; writing—original draft preparation, L.I.-P.; writing—review and editing, L.I.-P. and P.C.-M.; visualization, L.I.-P. and P.C.-M.; supervision, L.I.-P. and P.C.-M. Both authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Details of Models Applied to the Series

Appendix A.1. BDI Series

Table A1. ARMA(9,0) Model for the BDI series.
Table A1. ARMA(9,0) Model for the BDI series.
CoefficientStandard DeviationT Statisticp-Value
C−8.71 × 10−50.0006−0.13650.8915
AR(1)0.96950.013373.04170.0000
AR(2)−0.17910.0185−9.68530.0000
AR(3)−0.05460.0186−2.93230.0034
AR(4)0.06890.01873.69540.0002
AR(5)−0.04590.0187−2.46010.0139
AR(6)−0.00120.0187−0.06470.9484
AR(7)0.03300.01861.77180.0765
AR(8)0.01140.01850.61550.5382
AR(9)0.02410.01331.81260.0699
Table A2. GARCH(1,2) model for the BDI series.
Table A2. GARCH(1,2) model for the BDI series.
CoefficientStandard DeviationT Statisticp-Value
C6.97 × 10−63.62 × 10−719.25510.0000
RESID(−1)2 0.60100.025823.26630.0000
GARCH(−1)0.16380.02576.38140.0000
GARCH(−2)0.24730.018113.63030.0000
The linear model has been previously presented.
Table A3. EGARCH (2,3) model for the BDI series.
Table A3. EGARCH (2,3) model for the BDI series.
CoefficientStandard DeviationT Statisticp-Value
C(1)−0.01630.0021−7.81490.0000
C(2)0.49590.016030.96420.0000
C(3)−0.47720.0156−30.52520.0000
C(4)−0.00630.0010−6.26970.0000
C(5)1.47760.0087170.15380.0000
C(6)−0.21870.0003−675.20900.0000
C(7)−0.25920.0084−30.85260.0000
The linear model has been previously presented.

Appendix A.2. BDI1 Series

Table A4. ARMA(5,0) model for the BDI1 series.
Table A4. ARMA(5,0) model for the BDI1 series.
CoefficientStandard DeviationT Statisticp-Value
C9.98 × 10−60.00040.02480.9802
AR(1)0.74160.018140.98410.0000
AR(2)0.07490.02253.32370.0009
AR(3)0.02320.02261.02770.3042
AR(4)0.01230.02250.57500.5654
AR(5)−0.04390.0181−2.42870.0152
Table A5. GARCH(1,1) model for the BDI1 series.
Table A5. GARCH(1,1) model for the BDI1 series.
CoefficientStandard DeviationT Statistic p-Value
C1.46 × 10−62.06 × 10−77.074270.0000
RESID(−1)20.24350.013917.57750.0000
GARCH(−1)0.72460.011165.55490.0000
The linear model has been previously presented.
Table A6. EGARCH(3,1) model for the BDI1 series.
Table A6. EGARCH(3,1) model for the BDI1 series.
Coefficient Standard DeviationT Statisticp-Value
C(1)−0.18370.0213−8.63960.0000
C(2)0.41000.023517.43400.0000
C(3)−0.28750.0324−8.86990.0000
C(4)0.00610.02480.24410.8072
C(5)−0.02230.0053−4.20080.0000
C(6)0.99190.0016602.14250.0000
The linear model has been previously presented.

Appendix A.3. BDI2

Table A7. ARMA(3,0) model for BDI2 series.
Table A7. ARMA(3,0) model for BDI2 series.
CoefficientStandard ErrorT Statisticp-Value
C0.00130.00091.44320.1492
AR(1)1.04850.028237.18060.0000
AR(2)−0.11290.0409−2.75890.0059
AR(3)−0.13670.0282−4.84220.0000
Table A8. GARCH(1,1) model for BDI2 series.
Table A8. GARCH(1,1) model for BDI2 series.
CoefficientStandard DeviationT Statisticp-Value
C6.1 × 10−68.67 × 10−76.97440.0000
RESID(−1)20.50650.044011.49830.0000
GARCH(−1)0.45870.028915.89570.0000
The linear model has been previously presented.
Table A9. EGARCH(1,1) model for BDI2 series.
Table A9. EGARCH(1,1) model for BDI2 series.
CoefficientStandard DeviationT Statisticp-Value
C(5)−1.84850.2025−9.13010.0000
C(6)0.55420.034616.01380.0000
C(7)−0.02210.0235−0.93970.3474
C(8)0.86310.018347.16520.0000
The linear model has been previously presented.

Appendix A.4. BDI3 Series

Table A10. ARMA(7,0) model for BDI3 series.
Table A10. ARMA(7,0) model for BDI3 series.
CoefficientStandard DeviationT Statisticp-Value
C−0.00150.0021−0.72930.4659
AR(1)0.99640.026737.2930.0000
AR(2)−0.24340.0378−6.44110.0000
AR(3)−0.04360.0383−1.13760.2555
AR(4)0.09090.03832.37510.0177
AR(5)−0.03940.0384−1.02630.3049
AR(6)−0.03940.0378−1.04080.2982
AR(7)0.09650.02673.60850.0003
Table A11. GARCH(1,1) model for BDI3 series.
Table A11. GARCH(1,1) model for BDI3 series.
Coefficient Standard DeviationT Statisticp-Value
C4.83 × 10−53.85 × 10−612.54130.0000
RESID(−1) 20.39620.029613.40440.0000
GARCH(−1)0.45450.030115.10620.0000
The linear model has been previously presented.
Table A12. EGARCH model (1,1) for BDI3 series.
Table A12. EGARCH model (1,1) for BDI3 series.
Coefficient Standard DeviationT Statisticp-Value
C(1)−3.18770.2500−12.75200.0000
C(2)0.58200.036515.94240.0000
C(3)0.04540.02411.88610.0593
C(4)0.68130.027025.32000.0000
The linear model has been previously presented.

Appendix B. Results of Nonlinearity Tests

Appendix B.1. Results for Keenan, Tsay and White Tests

Table A13. Nonlinearity tests (BDI series).
Table A13. Nonlinearity tests (BDI series).
BDI
TestReturnsARMA(9,0)GARCH(1,2)EGARCH(2,3)
Keenan21.3281 [9]0.0009 [0]1.8963 [1]0.4427 [0]
3.956 × 10−60.97660.16850.5058
Tsay5.354 [9]0.0296 [0]3.013 [1]1.4550 [0]
7.235 × 10−280.99320.08260.2249
White70.4337132.428815.234632.896
5.551 × 10−16<2.2 × 10−160.00057.583 × 10−8
The first row shows the statistic associated with the test and in brackets, if necessary, the dimension of the model chosen to carry out the corresponding test. The second row shows the p-value. In bold, those values with a p-value < 0.05.
Table A14. Nonlinearity tests (subsamples).
Table A14. Nonlinearity tests (subsamples).
BDI1
TestReturnsARMA(5,0)GARCH(1,1)EGARCH(3,1)
Keenan11.3067 [5]2.8761 [26]2.3620 [5]0.1159 [0]
0.00080.09010.12440.7336
Tsay5.8340 [5]2.5410 [26]0.7093 [5]0.1371 [0]
4.519 × 10−122.978 × 10−390.77750.9379
White137.762019.13418.31535.1852
<2.2 × 10−167 × 10−50.000110.07483
BDI2
TestReturnsARMA(3,0)GARCH(1,1)EGARCH(1,1)
Keenan6.2973 [3]0.3993 [0]1.05581.4827
0.01220.52760.30440.2236
Tsay7.4420 [3]0.08871.2910 [5]1.349 [5]
7.631 × 10−80.96630.19990.1653
White46.4004155.32753.89035.8127
8.4 × 10−11<2.2 × 10−160.14300.0547
BDI3
TestReturnsARMA(7,0)GARCH(1,1)EGARCH(1,1)
Keenan8.8210 [7]0.0006 [0]0.0158 [0]0.0011 [0]
0.00300.98060.90010.9736
Tsay2.4470 [7]0.1001 [0]0.1445 [0]0.2253
4.16 × 1050.95990.93320.8788
White20.973542.737614.34539.373
2.79 × 10−55.2 × 10−100.00080.0092
Note: The first row shows the statistic associated with the test and in brackets, if necessary, the dimension of the model chosen to carry out the corresponding test. The second row shows the p-value.

Appendix B.2. Results of BDS Test

Appendix B.2.1. BDI Series

Table A15. BDS test results for ARMA(9,0) series.
Table A15. BDS test results for ARMA(9,0) series.
ϵm
0.5σ1.5σ
234.554736.265634.499831.6277
0.00000.00000.00000.0000
344.813242.653638.680034.5903
0.00000.00000.00000.0000
455.402447.591141.028135.6805
0.00000.00000.00000.0000
569.247852.753643.285436.6743
0.00000.00000.00000.0000
688.956058.584845.384737.4547
0.00000.00000.00000.0000
7116.322765.442647.597538.2830
0.00000.00000.00000.0000
8155.943573.690349.955939.1370
0.00000.00000.00000.0000
Note: The first row shows the statistic associated with the test. The second row shows the p-value. Values with a p-value less than the significance level 0.05 are shown in bold. m is the embedding dimension and ϵ is the distance expressed by the number of standard deviations of the data: (0.5−2)*σ.
Table A16. BDS test results for GARCH(1,2) series.
Table A16. BDS test results for GARCH(1,2) series.
ϵm
0.5σσ1.5σ
20.39350.64910.2717−0.4000
0.69390.51630.78580.6891
32.75722.67681.95760.9089
0.00580.00740.05030.3634
43.32523.18522.24570.9689
0.00090.00140.02470.3326
54.50274.11902.86311.3183
0.00000.00000.00420.1874
65.70094.98183.25561.4138
0.00000.00000.00110.1574
77.18005.86743.72911.6037
0.00000.00000.00020.1088
89.09226.84784.27201.8624
0.00000.00000.00000.0625
Note: The first row shows the statistic associated with the test. The second row shows the p-value. Values with a p-value less than the significance level 0.05 are shown in bold. m is the embedding dimension and ϵ is the distance expressed by the number of standard deviations of the data: (0.5–2)*σ.
Table A17. BDS test results for EGARCH(2,3) series.
Table A17. BDS test results for EGARCH(2,3) series.
ϵm
0.5σσ1.5σ
22.54573.09593.21822.6562
0.01090.00200.00130.0079
33.04703.68833.77213.1753
0.00230.00020.00020.0015
42.61073.43513.45302.9061
0.00900.00060.00060.0037
52.71393.43673.36062.8156
0.00660.00060.00080.0049
62.84673.44013.23642.6656
0.00440.00060.00120.0077
72.60443.34203.07512.4991
0.00920.00080.00210.0125
83.19093.26432.90282.3941
0.00140.00110.00370.0167
Note: The first row shows the statistic associated with the test. The second row shows the p-value. Values with a p-value less than the significance level 0.05 are shown in bold. m is the embedding dimension and ϵ is the distance expressed by the number of standard deviations of the data: (0.5–2)*σ.

Appendix B.2.2. BDI1 Series

Table A18. BDS test results for ARMA(5,0) series.
Table A18. BDS test results for ARMA(5,0) series.
ϵm
0.5σσ1.5σ
20.75860.89730.89760.5702
0.00000.00000.00000.0000
31.82681.81941.67491.1824
0.00000.00000.00000.0000
41.74201.92561.86921.3868
0.00000.00000.00000.0000
51.87582.24962.19121.6989
0.00000.00000.00000.0000
62.04732.59562.49541.9694
0.00000.00000.00000.0000
71.93522.72042.58352.0280
0.00000.00000.00000.0000
81.68522.84172.68662.1696
0.00000.00000.00000.0000
Note: The first row shows the statistic associated with the test. The second row shows the p-value. Values with a p-value less than the significance level 0.05 are shown in bold. m is the embedding dimension and ϵ is the distance expressed by the number of standard deviations of the data: (0.5−2)*σ.
Table A19. BDS test results for GARCH(1,1) series.
Table A19. BDS test results for GARCH(1,1) series.
ϵm
0.5σσ1.5σ
22.37601.75260.98850.4798
0.01750.07970.32290.6314
31.84351.0115−0.0020−0.6092
0.06530.31180.99840.5424
40.94050.1396−0.8142−1.3817
0.34700.88900.41550.1671
50.6708−0.2501−1.1556−1.6043
0.50230.80250.24790.1086
60.5129−0.2912−1.1979−1.6699
0.60800.77090.23100.0949
70.6479−0.3585−1.3044−1.7888
0.51700.72000.19210.0737
80.8758−0.3120−1.2992−1.7707
0.38110.75500.19390.0766
Note: The first row shows the statistic associated with the test. The second row shows the p-value. Values with a p-value less than the significance level 0.05 are shown in bold. m is the embedding dimension and ϵ is the distance expressed by the number of standard deviations of the data: (0.5–2)*σ.
Table A20. BDS test results for EGARCH(3,1) series.
Table A20. BDS test results for EGARCH(3,1) series.
ϵm
0.5σσ1.5σ
20.75860.89730.89760.5702
0.44810.36960.36940.5685
31.82681.81941.67491.1824
0.06770.06890.09400.2370
41.74201.92561.86921.3868
0.08150.05420.06160.1655
51.87582.24962.19121.6989
0.06070.02450.02840.0893
62.04732.59562.49541.9694
0.04060.00940.01260.0489
71.93522.72042.58352.0280
0.05300.00650.00980.0426
81.68522.84172.68662.1696
0.09190.00450.00720.0300
Note: The first row shows the statistic associated with the test. The second row shows the p-value. Values with a p-value less than the significance level 0.05 are shown in bold. m is the embedding dimension and ϵ is the distance expressed by the number of standard deviations of the data: (0.5−2)*σ.

Appendix B.2.3. BDI2 Series

Table A21. BDS test results for ARMA(3,0) series.
Table A21. BDS test results for ARMA(3,0) series.
ϵm
0.5σσ1.5σ
29.741211.148511.611812.0367
0.00000.00000.00000.0000
313.001313.613213.389013.2809
0.00000.00000.00000.0000
415.932115.025314.047613.5862
0.00000.00000.00000.0000
520.047716.584714.552413.7261
0.00000.00000.00000.0000
625.605718.355614.843613.6048
0.00000.00000.00000.0000
733.534620.700915.235113.4274
0.00000.00000.00000.0000
843.310423.292415.525813.1658
0.00000.00000.00000.0000
Note: The first row shows the statistic associated with the test. The second row shows the p-value. Values with a p-value less than the significance level 0.05 are shown in bold. m is the embedding dimension and ϵ is the distance expressed by the number of standard deviations of the data: (0.5−2)*σ.
Table A22. BDS test results for GARCH(1,1) series.
Table A22. BDS test results for GARCH(1,1) series.
ϵm
0.5σσ1.5σ
2−1.2666−1.7295−1.3671−0.9377
0.20530.08370.17160.3484
3−0.9555−1.5514−1.3344−1.0348
0.33930.12080.18210.3008
4−0.6801−1.1087−1.2029−1.0970
0.49640.26760.22900.2726
5−0.1696−0.6867−0.9729−0.9113
0.86530.49220.33060.3621
60.7995−0.2071−0.8367−0.9386
0.42400.83590.40280.3480
72.09720.3694−0.4836−0.7133
0.03600.71180.62870.4757
83.10210.6057−0.3686−0.6972
0.00190.54470.71240.4857
Note: The first row shows the statistic associated with the test. The second row shows the p-value. m is the embedding dimension and ϵ is the distance expressed by the number of standard deviations of the data: (0.5−2)*σ.
Table A23. BDS test results for EGARCH(1,1) series.
Table A23. BDS test results for EGARCH(1,1) series.
ϵm
0.5σσ1.5σ
2−0.2853−0.04120.36000.7764
0.77540.96710.71880.4375
3−0.1681−0.07200.28170.6315
0.86650.94260.77810.5277
4−0.2056−0.10420.12410.3303
0.83710.91700.90130.7412
5−0.1537−0.1343−0.01690.2078
0.87780.89320.98650.8354
6−0.08830.0099−0.1771−0.0278
0.92960.99210.85940.9778
70.71190.1617−0.1412−0.0612
0.47650.87150.88770.9512
80.49970.0522−0.3087−0.2840
0.61730.95830.75760.7764
Note: The first row shows the statistic associated with the test. The second row shows the p-value. m is the embedding dimension and ϵ is the distance expressed by the number of standard deviations of the data: (0.5−2)*σ.

Appendix B.2.4. BDI3 Series

Table A24. BDS test results for ARMA(7,0) series.
Table A24. BDS test results for ARMA(7,0) series.
ϵm
0.5σσ1.5σ
213.954412.828211.43009.3409
0.00000.00000.00000.0000
317.624515.041713.042610.7543
0.00000.00000.00000.0000
420.873116.253313.689511.4360
0.00000.00000.00000.0000
525.452417.736414.097511.5748
0.00000.00000.00000.0000
631.708619.492214.533311.5911
0.00000.00000.00000.0000
739.939121.364415.032211.6788
0.00000.00000.00000.0000
853.259223.523515.579911.7950
0.00000.00000.00000.0000
Note: The first row shows the statistic associated with the test. The second row shows the p-value. Values with a p-value less than the significance level 0.05 are shown in bold. m is the embedding dimension and ϵ is the distance expressed by the number of standard deviations of the data: (0.5−2)*σ.
Table A25. BDS test results for GARCH(1,1) series.
Table A25. BDS test results for GARCH(1,1) series.
ϵm
0.5σσ1.5σ
24.00303.43702.38171.0449
0.00010.00060.01720.2961
35.04593.63732.20010.6671
0.00000.00030.02780.5047
45.84263.73812.06310.4807
0.00000.00020.03910.6307
56.81343.96532.06930.3062
0.00000.00010.03850.7595
67.82714.26092.04090.1319
0.00000.00000.04130.8950
78.68324.51072.08550.1126
0.00000.00000.03700.9104
810.49404.85702.12550.1413
0.00000.00000.03350.8876
Note: The first row shows the statistic associated with the test. The second row shows the p-value. Values with a p-value less than the significance level 0.05 are shown in bold. m is the embedding dimension and ϵ is the distance expressed by the number of standard deviations of the data: (0.5−2)*σ.
Table A26. BDS test results for EGARCH(1,1) series.
Table A26. BDS test results for EGARCH(1,1) series.
ϵm
0.5σσ1.5σ
23.65843.32612.23370.8868
0.00030.00090.02550.3752
34.22273.50752.31680.8348
0.00000.00050.02050.4038
44.68013.66792.41751.0168
0.00000.00020.01560.3092
55.63384.00642.55851.0689
0.00000.00010.01050.2851
66.82574.39792.65161.0627
0.00000.00000.00800.2879
77.41754.68022.78161.1140
0.00000.00000.00540.2653
88.80025.00992.89321.2049
0.00000.00000.00380.2282
Note: The first row shows the statistic associated with the test. The second row shows the p-value. Values with a p-value less than the significance level 0.05 are shown in bold. m is the embedding dimension and ϵ is the distance expressed by the number of standard deviations of the data: (0.5−2)*σ.

Appendix C. Chaotic Component Study

Appendix C.1. Correlation Dimension Results

Table A27. Correlation dimension estimates.
Table A27. Correlation dimension estimates.
BDI
mReturns 1ARMA(9,0) 5GARCH(1,2) 3EGARCH(2,3) 1
10.860.740.910.96
21.311.381.811.92
31.671.962.702.88
42.002.493.593.85
52.313.004.454.82
62.883.475.315.79
73.153.936.127.06
83.424.376.958.95
BDI1
mReturns 12ARMA(5,0) 5GARCH(1,1) 2EGARCH(3,1) 3
10.970.920.950.97
21.931.821.911.93
32.882.722.892.91
43.813.593.863.83
54.784.434.834.76
65.395.215.835.71
75.656.036.686.43
85.586.836.938.21
BDI2
mReturns 6ARMA(3,0) 3GARCH(1,1) 2EGARCH(1,1) 1
10.990.850.970.98
21.981.661.931.97
33.032.42.842.93
44.053.093.724.13
55.453.714.84.99
66.044.245.695.63
74.484.665.584.70
84.325.198.726.81
BDI3
mReturns 5ARMA(7,0) 4GARCH(1,1) 2EGARCH(1,1) 1
10.940.970.970.97
21.861.931.931.94
32.762.892.902.89
43.733.873.843.85
54.734.784.664.89
65.775.765.895.68
76.837.687.586.44
87.127.476.906.46
Next to each applied model, the optimal value of the delay time according to the criterion of minimum information is indicated as a superscript; m is the embedding dimension.

Appendix C.2. Lyapunov Test Results

Table A28. Estimates of the Lyapunov exponents.
Table A28. Estimates of the Lyapunov exponents.
(L,m,q)Lyapunov Exponentsp-Value *Hypothesis
BDI
Returns(1,6,1)−0.17012.3 × 10−183H1
ARMA(9,0)(1,6,1)−0.49541.2 × 10−46H1
GARCH(1,2)(5,6,3)−0.61630.000H1
EGARCH(2,3)(3,6,4)−0.69021.6 × 10−28H1
BDI1
Returns(1,3,2)−0.17250.0000H1
ARMA(5,0)(2,6,5)−0.41110.0000H1
GARCH(1,1)(5,6,4)−0.49687.6 × 10−108H1
EGARCH(3,1)(5,6,4)−0.57056.8 × 10−128H1
BDI2
Returns(1,1,1)−0.14153.3 × 10−12H1
ARMA(3,0)(1,6,2)−0.40722.9 × 10−203H1
GARCH(1,1)(4,6,5)−0.46431.2 × 10−26H1
EGARCH(1,1)(1,6,5)−0.40128.7 × 10−98H1
BDI3
Returns(1,6,1)−0.16752.1 × 10−33H1
ARMA(7,0)(1,6,4)−0.44172.2 × 10−62H1
GARCH(1,1)(3,6,5)−0,47636.0 × 10−147H1
EGARCH (1,1)(3,5,5)−0.49701.1 × 10−54H1
Notes: * At 5% significance level, the null hypothesis of the existence of a chaotic component is rejected for those p-values less than 0.05. L is the time delay, m is the embedding dimension and q is the number of hidden units in the single-layer feedforward neural network.

Appendix C.3. MGRM Test Results

Table A29. MGRM test results.
Table A29. MGRM test results.
SeriesExponent (α)Standard Deviation
BDI
GARCH(1,2)0.27300.0094
EGARCH(2,3)0.30340.0123
BDI1
GARCH(1,1)0.22600.0054
EGARCH(3,1)0.25370.0077
BDI2
GARCH(1,1)0.24560.0091
EGARCH(1,3)0.22900.0089
BDI3
GARCH(1,1)0.23000.0102
EGARCH(1,1)0.24500.0089

References

  1. Faggini, M. Chaotic time series analysis in economics: Balance and perspectives. Chaos 2014, 24, 042101. [Google Scholar] [CrossRef] [Green Version]
  2. Lahmiri, S. A study on chaos in crude oil markets before and after 2008 international financial crisis. Phys. A 2017, 466, 389–395. [Google Scholar] [CrossRef]
  3. Guegan, D. Chaos in Economics and Finance. Annu. Rev. Control 2009, 33, 89–93. [Google Scholar] [CrossRef] [Green Version]
  4. Sandubete, J.E.; Escot, L. Chaotic signals inside some tick-by-tick financial time series. Chaos Solit. Fractals 2020, 137, 109852. [Google Scholar] [CrossRef]
  5. Gilmore, C.G. An examination of nonlinear dependence in exchange rates, using recent methods from chaos theory. Glob. Finance J. 2001, 12, 139–151. [Google Scholar] [CrossRef]
  6. Coto-Millan, P.; Inglada-Pérez, L. Testing for nonlinearity and chaos in liquid bulk shipping. Transp. Res. Proc. 2020, 48, 1605–1614. [Google Scholar] [CrossRef]
  7. UNCTAD. Review of Maritime Transport; United Nations Conference on Trade and Development: Geneva, Switzerland, 2013. [Google Scholar]
  8. Lin, A.J.; Chang, H.Y.; Hsiao, J.L. Does the Baltic Dry Index drive volatility spillovers in the commodities, currency, or stock markets? Transp. Res. E-Log. 2019, 127, 265–283. [Google Scholar] [CrossRef]
  9. Rodrigue, J.P.; Comtois, C.; Slack, B. The Geography of Transport Systems; Routledge: New York, NY, USA, 2013. [Google Scholar]
  10. Lin, F.; Sim, N.C.S. Trade, income and the Baltic Dry Index. Eur. Econ. Rev. 2013, 59, 1–18. [Google Scholar] [CrossRef] [Green Version]
  11. Guan, F.; Peng, Z.; Wang, K.; Song, X.; Gao, J. Multi-Step Hybrid Prediction Model of Baltic Supermax Index Based on Support Vector Machine. Neural Netw. World 2016, 26, 219–232. [Google Scholar] [CrossRef] [Green Version]
  12. Chen, S.; Meersman, H.; Van de Voorde, E. Forecasting spot rates at main routes in the dry bulk market. Marit. Econ. Logist. 2012, 14, 498–537. [Google Scholar] [CrossRef]
  13. Stopford, M. Maritime Economics, 3rd ed.; Routledge: London, UK, 2009. [Google Scholar]
  14. Han, Q.B.; Yan, G.; Ning, G.; Yu, B. Forecasting Dry Bulk Freight Index with Improved SVM. Math. Probl. Eng. 2014, 2014, 460684. [Google Scholar] [CrossRef]
  15. Zeng, Q.; Qu, C.; Ng, A.; Zhao, X. A new approach for Baltic Dry Index forecasting based on empirical mode decomposition and neural networks. Marit. Econ. Logist. 2016, 18, 192–210. [Google Scholar] [CrossRef]
  16. Geman, H.; Smith, W.O. Shipping markets and freight rates: An analysis of the Baltic Dry Index. J. Altern. Invest. 2012, 5, 98–109. [Google Scholar] [CrossRef] [Green Version]
  17. Goulielmos, A.M.; Psifia, M.E. Forecasting weekly freight rates for one-year time charter 65000 dwt bulk carrier, 1989–2008, using nonlinear methods. Marit. Policy Manag. 2009, 36, 411–436. [Google Scholar] [CrossRef]
  18. Goulielmos, A.M.; Psifia, M. A study of trip and time charter freight rate indices: 1968-2003. Marit. Policy Manag. 2007, 34, 55–67. [Google Scholar] [CrossRef]
  19. Goulielmos, A.M.; Giziakis, C.; Georgantzi, A. An application of nonlinear methods to the prediction of future freight rates, 2006–2008. Int. J. Ship. Trans. Log. 2012, 4, 78–106. [Google Scholar]
  20. McKenzie, M.D. Chaotic behavior in national stock market indices: New evidence from the close return test. Glob. Finance J. 2001, 12, 35–53. [Google Scholar] [CrossRef]
  21. BenSaïda, A.; Litimi, H. High level chaos in the exchange and index markets. Chaos Solitons Fract. 2013, 54, 90–95. [Google Scholar] [CrossRef]
  22. Su, X.; Wang, Y.; Duan, S.; Ma, J. Detecting chaos from agricultural product price time series. Entropy 2014, 16, 6415–6433. [Google Scholar] [CrossRef]
  23. Lahmiri, S. On fractality and chaos in Moroccan family business stock returns and volatility. Physica A 2017, 473, 29–39. [Google Scholar] [CrossRef]
  24. Brooks, C. Testing for nonlinearity in daily sterling exchange rates. Appl. Financ. Econ. 1996, 6, 307–317. [Google Scholar] [CrossRef]
  25. LeBaron, B. Chaos and Nonlinear Forecastability in Economics and Finance. Philos Trans A Math Phys Eng Sci 1994, 348, 397–404. [Google Scholar]
  26. Hsieh, D.A. Chaos and nonlinear dynamics: Application to financial markets. J. Financ. 1991, 46, 1839–1877. [Google Scholar] [CrossRef]
  27. Adrangi, B.; Chatrath, A.; Dhanda, K.K.; Raffieeb, K. Chaos in oil prices? Evidence from futures markets. Energy Econ. 2001, 23, 405–425. [Google Scholar] [CrossRef] [Green Version]
  28. Serletis, A. Is There Chaos in Economic Time Series? Can J. Econ. 1996, 29, 210–212. [Google Scholar] [CrossRef]
  29. Yousefpoor, P.; Esfahani, M.S.; Nojumi, H. Looking for systematic approach to select chaos tests. Appl. Math. Comput. 2008, 198, 73–91. [Google Scholar] [CrossRef]
  30. Caraiani, P. Testing for nonlinearity and chaos in economic time series with noise titration. Econ. Lett. 2013, 120, 192–194. [Google Scholar]
  31. Gunay, S.; Kaşkaloğlu, K. Seeking a Chaotic Order in the Cryptocurrency Market. Math. Comput. Appl. 2019, 24, 36. [Google Scholar] [CrossRef] [Green Version]
  32. Inglada-Perez, L. A Comprehensive Framework for Uncovering Nonlinearity and Chaos in Financial Markets: Empirical Evidence for Four Major Stock Market Indices. Entropy 2020, 22, 1435. [Google Scholar] [CrossRef]
  33. Santra, S.S.; Dassios, I.; Ghosh, T. On the Asymptotic Behavior of a Class of Second-Order Nonlinear Neutral Differential Equations with Multiple Delays. Axioms 2020, 9, 134. [Google Scholar] [CrossRef]
  34. Dassios, I.; Zimbidis, A.; Kontzalis, C. The delay effect in classical Samuelson’s model. J. Econ. Struct. 2014, 3, 7. [Google Scholar] [CrossRef] [Green Version]
  35. Zhu, P.; Xiang, Q. Oscillation criteria for a class of fractional delay differential equations. Adv. Differ. Equ. 2018, 2018, 403. [Google Scholar] [CrossRef]
  36. Agarwal, R.P.; Karakoç, F. A survey on oscillation of impulsive delay differential equations. Comput. Math. Appl. 2010, 60, 1648–1685. [Google Scholar] [CrossRef] [Green Version]
  37. Glass, D.S.; Jin, X.; Riedel-Kruse, I.H. Nonlinear delay differential equations and their application to modeling biological network motifs. Nat. Commun. 2021, 12, 1788. [Google Scholar] [CrossRef]
  38. Matilla-Garcia, M.; Ruiz Marin, M. A new test for chaos and determinism based on symbolic dynamics. J. Econ. Behav. Organ. 2010, 76, 600–614. [Google Scholar] [CrossRef] [Green Version]
  39. Box, G.E.P.; Jenkins, G.M. Time series analysis: Forecasting and control; Holden-Day: San Francisco, CA, USA, 1970. [Google Scholar]
  40. Engle, R.F.; Bollerslev, T. Modelling the persistence of conditional variances. Economet. Rev. 1986, 5, 1–50. [Google Scholar] [CrossRef]
  41. Nelson, D.B. Condicional Heteroskedasticity in Asset Returns: A New Aprproach. Econometrica 1991, 59, 347–370. [Google Scholar] [CrossRef]
  42. Inglada-Perez, L. Uncovering nonlinear dynamics in air transport demand. Int. J. Transp. Econ. 2016, 43, 33–66. [Google Scholar]
  43. Dickey, D.; Fuller, W. Distribution of the estimators for autoregressive time series with a unit root. J. Am. Stat. Assoc. 1979, 74, 426–431. [Google Scholar]
  44. Phillips, P.C.B.; Perron, P. Testing for Unit Roots in Time Series Regression. Biometrika 1988, 75, 335–346. [Google Scholar] [CrossRef]
  45. Kwiatkowski, D.; Phillips, P.C.B.; Schmidt, P.; Shin, Y. Testing the Null Hypothesis of Stationarity against the Alternative of a Unit Root. J. Econom. 1992, 54, 159–178. [Google Scholar] [CrossRef]
  46. Keenan, D.M.A. Tukey nonadditivity-type test for time series nonlinearity. Biometrika 1985, 72, 39–44. [Google Scholar] [CrossRef]
  47. Tsay, R.S. Nonlinearity tests for time series. Biometrika 1986, 73, 461–466. [Google Scholar] [CrossRef]
  48. White, H. An additional hidden unit test for neglected nonlinearity in multilayer feedforward networks. In Proceedings of the International Joint Conference on Neural Networks, New York, NY, USA, 1–3 July 1989; IEEE Press: New York, NY, USA, 1989; pp. 451–455. [Google Scholar]
  49. Brock, W.A.; Dechert, W.D.; Sheinkman, J.A. A Test of Independence Based on the Correlation Dimension; SSRI no. 8702; Department of Economics, University of Wisconsin: Madison, WI, USA, 1987. [Google Scholar]
  50. Brock, W.A.; Scheinkman, J.A.; Dechert, W.D.; LeBaron, B. A test for independence based on the correlation dimension. Econom. Rev. 1996, 15, 197–235. [Google Scholar] [CrossRef]
  51. Grassberger, P.; Procaccia, I. Dimensions and entropies of strange attractors from a fluctuating dynamics approach. Physica 1984, 13, 34–54. [Google Scholar] [CrossRef]
  52. Eckmann, J.; Ruelle, D. Ergodic theory of chaos and strange attractors. Rev. Mod. Phys. 1985, 57, 617–656. [Google Scholar] [CrossRef]
  53. Schwarz, G. Estimating the dimension of a model. Ann. Stat. 1978, 6, 461–464. [Google Scholar] [CrossRef]
  54. Barkoulas, J.T.; Chakraborty, A.; Ouandlous, A. A metric and topological analysis of determinism in the crude oil spot market. Energy Econ. 2012, 34, 584–591. [Google Scholar] [CrossRef]
  55. Fraser, A.M.; Swinney, H.L. Independent Coordinates for Strange Attractors from Mutual Information. Phys. Rev. A 1986, 33, 1134–1140. [Google Scholar] [CrossRef]
  56. Takens, F. Detecting strange attractors in turbulence. In Dynamical Systems and Turbulence; Rand, D., Young, L., Eds.; Springer: Berlin/Heidelberg, Germany, 1981; pp. 366–381. [Google Scholar]
  57. Abraham, C.; Biau, G.; Cadre, B. Chaotic properties of mappings on a probability space. J. Math. Anal. Appl. 2002, 266, 420–431. [Google Scholar] [CrossRef] [Green Version]
  58. Gu, R.; Chen, X.; Li, X. Chaos recognition and fractal analysis in the term structure of Shanghai Interbank Offered Rate. Phys. A 2014, 412, 101–112. [Google Scholar] [CrossRef]
  59. Barrio, R.; Martínez, M.A.; Serrano, S.; Wilczak, D. When chaos meets hyperchaos: 4D Rössler model. Phys. Lett. A 2015, 379, 2300–2305. [Google Scholar] [CrossRef]
  60. Nychka, D.; Ellner, S.; Gallant, A.R.; McCaffrey, D. Finding chaos in noisy systems. J. R. Stat. Soc. Series B Stat. Methodol. 1992, 54, 399–426. [Google Scholar] [CrossRef]
  61. Bask, M.; Tung, L.; Widerberg, A. The stability of electricity prices: Estimation and inference of the Lyapunov exponents. Phys. A 2007, 376, 565–572. [Google Scholar] [CrossRef] [Green Version]
  62. Lahmiri, S. Investigating existence of chaos in short and long term dynamics of Moroccan exchange rates. Phys. A 2017, 465, 655–661. [Google Scholar] [CrossRef]
  63. Wolf, A.; Swift, J.B.; Swinney, H.L.; Vastano, J.A. Determining Lyapunov exponents from a time series. Phys. D 1985, 16, 285–317. [Google Scholar] [CrossRef] [Green Version]
  64. Shintani, M.; Linton, O. Nonparametric neural network estimation of Lyapunov exponents and direct test for chaos. J. Econom. 2004, 120, 1–33. [Google Scholar] [CrossRef] [Green Version]
  65. McCaffrey, D.; Ellner, S.; Gallant, A.; Nychka, D. Estimating the Lyapunov exponent of a chaotic system with nonparametric regression. J. Am. Stat. Assoc. 1992, 87, 682–695. [Google Scholar] [CrossRef]
  66. Bailey, B.A.; Ellner, S.; Nychka, D.W. Chaos with confidence: Asymptotics and applications of local Lyapunov exponents. In Nonlinear Dynamics and Time Series: Building a Bridge between the Natural and Statistical Sciences; Cutler, C.D., Kaplan, D.T., Eds.; Fields Institute Communications, American Mathematical Society: Providence, RI, USA, 1997; pp. 115–133. [Google Scholar]
  67. Eckmann, J.P.; Kamphorst, S.O.; Ruelle, D. Recurrence plots of dynamical systems. Europhys. Lett. 1987, 4, 973–977. [Google Scholar] [CrossRef] [Green Version]
  68. Cajueiro, D.O.; Tabak, B.M. Long-range dependence and multifractality in the term structure of LIBOR interest rates. Phys. A 2007, 73, 603–614. [Google Scholar] [CrossRef]
  69. Stavroyiannis, S.; Babalos, V.; Bekiros, S.; Lahmiri, S.; Uddin, G.S. The high frequency multifractal properties of Bitcoin. Phys. A 2019, 520, 62–71. [Google Scholar] [CrossRef]
  70. Fama, E.F. Efficient Capital Markets: A Review of Theory and Empirical Work. J. Financ. 1970, 25, 383–417. [Google Scholar] [CrossRef]
  71. Jing, L.; Marlow, P.B.; Hui, W. An analysis of freight rate volatility in dry bulk shipping markets. Marit. Policy Manag. 2008, 35, 237–251. [Google Scholar] [CrossRef]
  72. MacKinnon, J.G. Numerical distribution functions for unit root and cointegration tests. J. Appl. Econom. 1996, 11, 601–618. [Google Scholar] [CrossRef] [Green Version]
  73. Barnett, W.A.; Gallant, A.R.; Hinich, M.J.; Jungeilges, J.A.; Kaplan, D.T.; Jensen, M.J. A single-blind controlled competition among tests for nonlinearity and chaos. J. Econom. 1997, 82, 157–192. [Google Scholar] [CrossRef] [Green Version]
  74. Kennel, M.B.; Brown, R.; Mand Abarbanel, H.D. Determining embedding dimension for phase-space reconstruction using a geometrical construction. Phys. Rev. A 1992, 45, 3403. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Plot of four Baltic Dry Index (BDI) series. Time points (year) are on x-axis and observations are on y-axis.
Figure 1. Plot of four Baltic Dry Index (BDI) series. Time points (year) are on x-axis and observations are on y-axis.
Mathematics 09 02065 g001
Figure 2. Evolution of BDI returns.
Figure 2. Evolution of BDI returns.
Mathematics 09 02065 g002
Figure 3. (A) Recurrence plot for the returns and ARMA series from the BDI series; (B) Recurrence plot for the GARCH and EGARCH series from the BDI series.
Figure 3. (A) Recurrence plot for the returns and ARMA series from the BDI series; (B) Recurrence plot for the GARCH and EGARCH series from the BDI series.
Mathematics 09 02065 g003aMathematics 09 02065 g003b
Figure 4. (A) Recurrence plot for the returns and ARMA series from the BDI1 series; (B) Recurrence plot for the GARCH and EGARCH series from the BDI1 series.
Figure 4. (A) Recurrence plot for the returns and ARMA series from the BDI1 series; (B) Recurrence plot for the GARCH and EGARCH series from the BDI1 series.
Mathematics 09 02065 g004aMathematics 09 02065 g004b
Figure 5. (A) Recurrence plot for the returns and ARMA series from the BDI2 series; (B) Recurrence plot for the GARCH and EGARCH series from the BDI2 series.
Figure 5. (A) Recurrence plot for the returns and ARMA series from the BDI2 series; (B) Recurrence plot for the GARCH and EGARCH series from the BDI2 series.
Mathematics 09 02065 g005aMathematics 09 02065 g005b
Figure 6. (A) Recurrence plot for the returns and ARMA series from the BDI3 series; (B) Recurrence plot for the GARCH and EGARCH series from the BDI3 series.
Figure 6. (A) Recurrence plot for the returns and ARMA series from the BDI3 series; (B) Recurrence plot for the GARCH and EGARCH series from the BDI3 series.
Mathematics 09 02065 g006
Table 1. Techniques used.
Table 1. Techniques used.
TechniqueReferenceFeature Tested
Augmented Dickey Fuller (ADF)[43]Stationarity (unit roots)
Phillips–Perron (PP)[44]Stationarity (unit roots)
Kwiatkowski–Phillips–Schmidt–Shin (KPSS)[45]Stationarity (unit roots)
Keenan[46]Nonlinearity
Tsay[47]Nonlinearity
White[48]Nonlinearity
Brock, Dechert and Sheinkman test-BDS[49,50]Independence (iid); Randomness; Nonlinearity
Correlation Dimension[51]Chaos
Lyapunov Exponent[21]Chaos
MGRM[38]Chaos
Recurrence Plots[52]Chaos
Note: MGRM stands for the test proposed by Matilla-Garcia and Ruiz Marin [38]. Source: [32] and own work.
Table 2. Descriptive statistics.
Table 2. Descriptive statistics.
BDI Series
Original SeriesReturnsARMA(9,0)GARCH(1,2)EGARCH(2,3)
Mean2283.77−7.2 × 10−54.26 × 10−110.01500.0078
Standard Deviation19090.01480.00840.97231.0030
Median1531−9.3 × 10−5−0.0001−0.0191−0.0282
Minimum647−0.1200−0.0783−5.8461−5.4268
Maximum11,7930.14000.097812.61057.1285
Skewness Coefficient2.41−0.06000.37661.11980.4531
Kurtosis Coefficient9.2112.800018.572815.13646.8089
Jarque–Bera Test14,634.4 *22,759.8 *57,559.1 *36,065.6 *3629.8 *
N56925691568356835683
BDI1 Series
Original SeriesReturnsARMA(5,0)GARCH(1,1)EGARCH(3,1)
Mean1360.462.7 × 10−5−3 × 10−11−0.00510.0147
Standard Deviation314.980.00720.00420.96341.0022
Median1353−0.0006−0.0001−0.0399−0.0152
Minimum776.00−0.03900.0385−4.9043−4.7600
Maximum23520.0430−0.02808.58577.1635
Skewness Coefficient0.49060.58700.73820.62850.4990
Kurtosis Coefficient3.10995.901911.67368.10906.8966
Jarque–Bera Test124.23 *1248.6 *9847.4 *3521.2 *2058.2 *
N30593058305330533053
BDI2 Series
Original SeriesReturnsARMA(3,0)GARCH(1,1)EGARCH(1,1)
Mean4157.690.00133.2 × 10−12−0.0490−0.0235
Standard Deviation2028.520.01340.00660.96110.9997
Median4038.500.0017−0.0002−0.0885−0.0668
Minimum1530−0.1200−0.0444−5.3013−4.9095
Maximum11,793−0.14000.08554.27584.1532
Skewness Coefficient1.44470.28191.52000.01660.0341
Kurtosis Coefficient5.24474.300029.83874.94904.6545
Jarque–Bera Test690.58 *100.91 *37,632196.0141.44
N12381238123812381238
BDI3 Series
Original SeriesReturnsARMA(7,0)GARCH(1,1)EGARCH(1,1)
Mean2645.20−0.0015−8.90 × 10−110.00360.0066
Standard Deviation2427.370.02510.01420.96491.0004
Median1817.50−0.0019−0.0006−0.0603−0.0600
Minimum647−0.1200−0.0767−4.3534−4.5133
Maximum11,7930.13700.09865.80116.2652
Skewness Coefficient1.89800.01200.23600.68570.6806
Kurtosis Coefficient5.94156.08407.71697.14887.2499
Jarque–Bera Test1341.42 *553.2 *1307.11110.581158.3
N13961396139613961396
Note: Statistics for the considered time series: original series, return series, residuals of autoregressive moving average (ARMA) process, as well as standardized residuals from generalized autoregressive conditional heteroscedastic (GARCH) model and exponential GARCH (EGARCH) process. Skewness corresponds to Fisher’s skewness coefficient. *: Asterisk denotes significance at the 5% level. N is the number of observations. Source: Own work.
Table 3. Stationarity analysis.
Table 3. Stationarity analysis.
BDIOriginal SeriesReturn Series
1. Augmented Dickey–Fuller Test
Constant−3.3632 (0.0124)−18.5200 (0.0000)
Constant and Linear Trend−3.4880 (0.0407)−18.5215 (0.0000)
2. Phillips–Perron Test
Constant−2.7104 (0.0723)−23.4660 (0.0000)
Constant and Linear Trend−2.7660 (0.2101)−23.4661 (0.0000)
3. Kwiatkowski–Phillips–Schmidt–Shin Test
Constant2.07040.0510
Constant and Linear Trend0.58540.0385
BDI1Original SeriesReturn Series
1. Augmented Dickey–Fuller Test
Constant−2.7328 (0.0686)−16.0082 (0.0000)
Constant and Linear Trend−2.6569 (0.2549)−16.0196 (0.0000)
2. Phillips–Perron Test
Constant−2.2710 (0.1816)−18.3371 (0.0000)
Constant and Linear Trend−2.1487 (0.5176)−18.3492 (0.0000)
3. Kwiatkowski–Phillips–Schmidt–Shin Test
Constant0.89110.0892
Constant and Linear Trend0.25900.0593
BDI2Original SeriesReturn Series
1. Augmented Dickey–Fuller Test
Constant−0.0856 (0.9491)−13.1635 (0.0000)
Constant and Linear Trend−0.7875 (0.9653)−13.1591 (0.0000)
2. Phillips–Perron Test
Constant−0.2200 (0.9334)−8.5187 (0.0000)
Constant and Linear Trend−0.9141 (0.9528)−8.5131 (0.0000)
3. Kwiatkowski–Phillips–Schmidt–Shin Test
Constant1.61400.1329
Constant and Linear Trend0.59250.1274
BDI3Original SeriesReturn series
1. Augmented Dickey–Fuller Test
Constant−2.5787 (0.0976)−14.5984 (0.0000)
Constant and Linear Trend−2.3326 (0.4154)−14.6130 (0.0000)
2. Phillips–Perron Test
Constant−2.6358 (0.0860)−12.2124 (0.0000)
Constant and Linear Trend−2.6044 (0.2784)−12.2111 (0.0000)
3. Kwiatkowski–Phillips–Schmidt–Shin Test
Constant2.34570.0690
Constant and Linear Trend0.30920.0385
Note: Augmented Dickey–Fuller p-values correspond to [72] one-sided p-values. For the Phillips–Perron test, lags were based on bandwidth Newey–West using the Bartlett kernel. Critical values for the Kwiatkowski, Phillips, Schmidt and Shin test are 0.463 and 0.146, respectively, for the constant and linear plus linear trend model. Significant values at the 95% confidence level are in bold.
Table 4. Classes of models applied to each series.
Table 4. Classes of models applied to each series.
SeriesLinear ModelNonlinear Model
[ARMA][GARCH][EGARCH]
BDIARMA(9,0)GARCH(1,2)EGARCH(2,3)
BDI1ARMA(5,0)GARCH(1,1)EGARCH(3,1)
BDI2ARMA(3,0)GARCH(1,1)EGARCH(1,1)
BDI3ARMA(7,0)GARCH(1,1)EGARCH(1,1)
Table 5. Summary of results for the nonlinear analysis.
Table 5. Summary of results for the nonlinear analysis.
BDI Series
TestHoReturnsARMA(9,0)GARCH(1,2)EGARCH(2,3)
KeenanLaRNRNRNR
TsayLaRNRNRNR
WhiteLaRRRR
BDSLb-RMixtureR
BDI1 Series
TestHoReturnsARMA(5,0)GARCH(1,1)EGARCH(3,1)
KeenanLaRNRNRNR
TsayLaRRNRNR
WhiteLaRRRNR
BDSLb-RNRNR
BDI2 Series
TestHoReturnsARMA(3,0)GARCH(1,1)EGARCH(1,1)
KeenanLaRNRNRNR
TsayLaRNRNRNR
WhiteLaRRNRR
BDSLb-RNRNR
BDI3 Series
TestHoReturnsARMA(7,0)GARCH(1,1)EGARCH(1,1)
KeenanLaRNRNRNR
TsayLaRNRNRNR
WhiteLaRRRR
BDSLb-RRR
Ho: Null hypothesis. La: Linear in mean. Lb: This test can be applied as a linear test once the linear dependence has been eliminated. R represents that null hypothesis is rejected; NR represents that null hypothesis is not rejected. Mixture represents that the test reported no clear results.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Inglada-Pérez, L.; Coto-Millán, P. A Chaos Analysis of the Dry Bulk Shipping Market. Mathematics 2021, 9, 2065. https://doi.org/10.3390/math9172065

AMA Style

Inglada-Pérez L, Coto-Millán P. A Chaos Analysis of the Dry Bulk Shipping Market. Mathematics. 2021; 9(17):2065. https://doi.org/10.3390/math9172065

Chicago/Turabian Style

Inglada-Pérez, Lucía, and Pablo Coto-Millán. 2021. "A Chaos Analysis of the Dry Bulk Shipping Market" Mathematics 9, no. 17: 2065. https://doi.org/10.3390/math9172065

APA Style

Inglada-Pérez, L., & Coto-Millán, P. (2021). A Chaos Analysis of the Dry Bulk Shipping Market. Mathematics, 9(17), 2065. https://doi.org/10.3390/math9172065

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop