ARS MATHEMATICA CONTEMPORANEA
Also available at http://amc-journal.eu ISSN 1855-3966 (printed edn.), ISSN 1855-3974 (electronic edn.)
ARS MATHEMATICA CONTEMPORANEA 13 (2017) 63-79
The JLS model with ARMA/GARCH errors
Faculty of Mathematics and Physics, University of Ljubljana, Jadranska 19, SI-1000 Ljubljana, Slovenia
Matjaž Omladic
Institute of Mathematics, Physics and Mechanics, Jadranska 19, SI-1000 Ljubljana, Slovenia
Received 16 October 2014, accepted 3 December 2014, published online 21 October 2016
Prior to crashes, the stock index price time series is characterised by the Log-Periodic Power Law (LPPL) equation of the Johansen-Ledoit-Sornette (JLS) model, which leads to a critical point indicating a change to a new market regime. In this paper, we describe the hierarchical diamond lattice, upon which the JLS model is derived, using the diamond lattice operation Di and derive the recursion for the coefficients of the growth function in a diamond lattice rooted at the main root vertex rm. Further, to verify the adequacy of the JLS model, we analyse the model's residuals and propose its generalization, using the ARMA/GARCH error model. We determine the ARMA/GARCH orders using the extended autocorrelation function (EACF) method and compare the results with those of the Akaike and Bayesian Information Criteria. Using the data for 33 major world stock indices we show, that proposed generalization of the JLS model in general performs better in predicting the market regime changes and has also the ability to recognise false bubble identification, indicated by the JLS model.
Keywords: Graph operations, hierarchical diamond lattice, JLS model, financial bubbles and crashes, ARMA/GARCH errors.
Math. Subj. Class.: 05C76, 05C10, 82B20, 62P20, 62M10, 91B84
*Xlab d.o.o., Pot za Brdom 100, SI-1000 Ljubljana, Slovenia
E-mail addresses: spela.jezernik@gmail.com (Spela Jezernik Sirca), matjaz@omladic.net (Matjaž Omladic)
Spela Jezernik Sirca

Abstract
©® This work is licensed under http://creativecommons.org/licenses/by/3.0/
64
Ars Math. Contemp. 13 (2017) 125-136
1 Introduction
Financial bubbles and crashes are fascinating events for academics and practitioners alike, and such occurrences are especially interesting in the field of econophysics. These events are extremely important because of their usually strong impacts on not only financial markets but also the global economy. Unfortunately, there is still no general agreement in the literature on what defines a financial bubble or crash. A financial bubble may be recognised as a long-term positive deviation of a financial asset's market price from its fundamental value. A crash may be defined as a sudden dramatic decline of market price over a short time period. Understanding the behaviour of financial markets and the relationship between financial bubbles and crashes may help to minimise the damage of the speculative bubbles that end up with crashes. Consequently, identifying a financial bubble and predicting its end has become very important issue in financial markets behaviour research.
The Johansen-Ledoit-Sornette (JLS) model [14] was developed to describe the dynamics of financial bubbles and crashes. The model states that, prior to crashes, the mean function of a stock index price time series is characterised by a Log-Periodic Power Law (LPPL) equation that leads to a critical point indicating the change to a new market regime, either a large crash or a change in the average growth rate. The model assumes the presence of two types of agents in the market, namely a group of agents with rational expectations and a group of irrational agents with herding behaviours, and these agents potentially lead to the development of speculative bubbles. These agents are organised into a network in which each exists in only one of two possible states (buy or sell), while their trading actions depend on the opinions of other agents and on external influences. If the tendency of irrational agents to imitate their neighbours increases up to a certain critical point, then a large proportion of agents will be in the same state (sell) at the same time, thus causing a crash. In the JLS model, bubbles are characterised by faster-than-exponential price growth due to herding behaviours and imitation of irrational agents during the bubble period. The key parameter of the model is the critical time tc, which is interpreted as the moment at which the bubble ends and the transition to another market regime begins.
Numerous empirical results have been reported by several authors on this subject. The JLS model has been used in various types of markets, such as the bubbles of stock market indices [14, 16, 9, 44, 12], the anti-bubbles in different financial markets [15, 32, 40], exchange rate bubbles [20], the oil bubble [31], real estate bubbles [42, 43], corporate bond spread bubbles [3], credit default swap (CDS) spread bubbles [36], and the repo market size [38]. Most of the published research papers on the JLS model have focused on the existence of log-periodic fluctuations by fitting LPPL equation to the data. Although some papers have included several statistical methods for the detection of log-periodicity [14,29,6, 39,41,2], only a few papers have focused on the JLS model residuals [8,19,12].
The aim of this paper is to propose the JLS model generalization, based on the analysis of the JLS model residuals. Specifically, we investigate the presence of ARMA/GARCH patterns in the JLS model residuals, wherein we also compare our results with the log-peri-odic-AR(1)-GARCH(1,1) specification, proposed in [8]. According to ARMA/GARCH model determination, we examine the adequacy of the JLS model. In doing so, we compare the critical time parameter estimates, calculated with the JLS model, versus those calculated with generalized JLS model. We explore whether the generalized JLS model improves the JLS model estimates of the critical time parameters. To assure the generality of the results, we perform an analysis on large number of data samples.
The rest of the paper is organised as follows. Section 2 describes derivation of the
S. J. Sirca and M. Omladic: The JLS model with ARMA/GARCH errors
65
JLS model, together with the q-state Potts model and growth functions on hierarchical diamond lattice. In this section we also present details of the JLS model estimation, optimisation and verification. Our proposed generalization of the JLS model, methodology for the ARMA/GARCH model determination and parameter estimation are described in section 3. The data, empirical results of our analysis and main contributions of this paper are reported in section 4. Section 5 concludes the paper.
2 The Johansen-Ledoit-Sornette model 2.1 Motivation
Financial markets consist of numerous interacting traders, that differ in size from small individuals to large institutional traders, such as pension funds. Moreover, all traders worldwide are organised inside a social network (family, friends, etc), within which they locally influence each other. The structure of financial markets resembles to hierarchical systems with traders on all different levels of the market.
To develop the Johansen-Ledoit-Sornette (JLS) model, a model of rational imitation, Johansen et al. [14] used hierarchical diamond lattice representation for the structure of financial markets. In the case of hierarchical diamond lattice discussed by Berker and Ostlund [1], the lattice is generated in an iterative manner as shown in Figure 1. This is quite realistic model of complicated network of interactions between traders. The model is derived using the q-state Potts model on hierarchical diamond lattice defined in [5], where free energy exhibits log-periodic oscillations as the critical point is approached. For more details, see [13].
IO
Figure 1: First few steps of building a hierarchical diamond lattice.
2.2 The q-state Potts model on hierarchical diamond lattice
Let G be a graph and consider a set {1,2,..., q} of q elements, called spins. A state of a graph G is an assignment of a single spin to each vertex of the graph. Denote by V (G) = {v\,..., vn} vertex set of G and by E (G) edge set of G. Then the state of G is a function a : V (G) ^ {1,2,..., q}. Let S (G) denote the set of states of G.
The interaction energy may be thought of simply as weights on the edges of the graph G. Denote by Je = J^j = JVi,Vj interaction energy on an edge e = {vj, Vj}. Then the Hamiltonian is
h (a) = -	Jj,jS (aj, aj),
{i,j}eE(G)
where a is a state of graph G, ai is the spin at vertex vi, S is Kronecker delta function and
66
Ars Math. Contemp. 13 (2017) 125-136
each edge {i, j} is assigned an interaction energy Jjj.
Let the spins be positioned on a hierarchical diamond lattice constructed by the iterative process shown in Figure 1. Denote by Gn graph on the nth level of hierarchical diamond lattice. Then the q-state Potts model partition function on the nth level is defined as
Zn (Gn)= £ e-^.
aeS(Gn)
The Potts model partition function is the sum of all possible states of an exponential function of the Hamiltonian.
There exists a connection between the Potts model, which is useful to study phase transitions and critical phenomena in physics, and the graph theory, for example the Tutte polynomial [35] or Chromatic polynomial [28]. The graph theory is mathematical area, useful to describe and study the relations between participants in networks, such as physical, biological, social or economic networks. For some basic concepts used in models of economic networks, see [17].
Figure 1 depicts the first few steps of the diamond lattice operation, denoted by Di, where the original graph is a single edge. The graphs in the Figure are G0 = K2, Gi = Di(G0), G2 = Di(Gi), G3 = Di(G2). The original graph is planar, so is each next graph in the sequence. Di can be applied on any map (graph embedded on a surface). Note that Di is a composite operation. In [7, 24] several operations on maps are considered. One, Pa, parallelization replaces each edge by a pair of parallel edges and another one Su1, one-dimensional subdivision subdivides each edge of the original map. In this way Di(M) = Su1(Pa(M)). Note that Di may be regarded as an operation on map or due to its simplicity also as an operation on the underlying graph. A theory of representations of graphs and maps has been laid down by Pisanski and Zitnik [26], where such operations were considered. Repeated operations were used in other contexts, see for instance [25]. Operations on maps have been studied in connection with symmetry in several papers, see for instance [23, 22, 4, 10].
2.3 Growth function in rooted diamond lattice
Recall the definition of growth function in rooted graphs [25]. Let G be a connected, finite or locally finite (infinite graph with finite vertex degrees) graph. Let V(G) be the set of vertices of G and let r e V(G) be the so-called root of graph G. The distance d(u,v) between u, v e V(G) is defined as number of edges on the shortest path between u and v. Further, we define the (spherical) growth sequence as {¿(G, r, n)|n = 0,1,2,... }, where S(G, r, n) denote the number of vertices at distance n from root r in graph G. Then the (spherical) growth function of graph G rooted at r can be written as follows:
f (G,r,x) = £ ¿(G,r,n)xn,
n=0
i.e. the generating function for S(G, r, n) of graph G at root r. Here we limit to the case when the root is single vertex, vertex root, although we can extend the definition by allowing that a root is any induced subgraph of graph G.
We start the construction of hierarchical diamond lattice with a graph G0 = K2. Its growth function is f (G0, r, x) = 1 + x and is independent of the selected root r. Using the diamond lattice operation on graph G0 we get graph Gi = Di(G0), for which
S. J. Sirca and M. Omladic: The JLS model with ARMA/GARCH errors
67
f (Gi, r, x) = 1 + 2x + x2 is also independent of r. If we use diamond lattice operation at least two times, we have to calculate each of the growth functions separately depending on the root. Here we first limit ourselves to the case of selecting one of the two vertices at the top and bottom in each graph Gn, called main root vertex, denoted by rm (white vertices of graph Gn on Figure 1). Note that since each diamond lattice Gn can be mirrored horizontally (vertically for n > 0), horizontally (vertically) symmetrically selected roots produce an identical growth function. In the next subsection we examine the cases of selecting vertex root different from rm in Gn for n = 1,2, 3,4. The results for all possible growth functions of a graph Gn for n > 4 will be published elsewhere.
Let us present the growth function of graph Gn rooted at rm for n = 2 and n = 3:
f (G2, rm, x) = 1 + 4x + 2x2 + 4x3 + x4
f (G3,rm,x) = 1 + 8x + 4x2 + 8x3 + 2x4 + 8x5 + 4x6 + 8x7 + x8.
If we denote the growth function of graph Gn rooted at rm by gn(x), we can write the system of equations that determine the growth functions f (Gn, rm, x) as follows:
go(x)	=	(1 + x)
gi(x)	=	(1+ x)(2go(x) - (1 + x))
g2(x)	=	(1+ x2)(2gi(x) - (1 + x2))
53(x)	=	(1+ x4)(2g2(x) - (1 + x4))
gn(x) = (1+ X2" )(2gn_i(x) - (1 + :

(2.1)
Moreover, if we define the function
/o(x) = 1 and /n(x) = ^Q(1 + x2") for n > 0,
i=i
(2.2)
we can write the recursion of the growth function /(Gn, rm, x) also as:
gn(x) = gn_l(x2) + 2"x/„_i(x).
To derive the recursion for the list of coefficients of gn (i.e. for the growth sequence of graph Gn rooted at rm), we first define some operations on lists. Let w = {wi, w2,..., wm} and w = {wW1, w2,..., wn} be two lists of coefficients. Let a(k) = {a, a,..., a} denote the list of k repetitions of the value a. Then we define the following operations on the lists:
t(w, w m(w, k v(w, a b(w c(w
= w + w = {wi, w2, . . . , wm, wi, w2, . . . , wn} = w • k = {wik, w2k,. .., wmk}
=wa
(m_i) =
= { wi , a, w2 ,
,wm_ i, a, w„
w[ 1] = {w2, w3, . . . , wm} w[ 1, m] = {w2, w3, .. ., wm_i} .
Let an denote the list of coefficients of gn. Then we derive the following rule for generating the list of coefficients for gn from the list of coefficients for gn_i and g0:
ao = {1, 1}
an = t(6(ao),t(m(t(6(a„_i),c(a„_i)), 2),b(ao))).
(2.3)
2
}
68
Ars Math. Contemp. 13 (2017) 125-136
Alternatively, the recursion in (2.3) can be written as follows:
a„ = v(a„_i, 2n).	(2.4)
The first four lists an of coefficients for gn are:
ao = {1, 1}
ai = {1,2,1}
a2 = {1, 4, 2,4,1}
a3 = {1, 8,4, 8, 2, 8, 4, 8,1}
At each step we get a list an of coefficients for gn of length 2n + 1.
2.4 Different growth functions of Gn for n = 1, 2, 3, 4
In the previous subsection we examined the growth function in diamond lattice Gn rooted at rm. Here we examine the number of different growth functions in Gn and calculate the growth functions separately depending on the root r different from rm in Gn for n = 1,2,3,4 (the graph G0 has only two main root vertices).
For a graph G and v e V(G), the degree of a vertex v is denoted by dG(v) and the degree set of a graph G (i.e. the set of (distinct) degrees of the vertices of G) by D(G), where we write the degree set in ascending order. Then, we define the set VD(G) = {pki(G)|kj e D(G), i = 1, 2,..., |D(G)|}, where pk(G) denotes the number of vertices of a degree k in G. Let A be the automorphism group of G that partitions the vertex set V(G) into orbits:
V(G) = [vi] U [v2] U • • • U [vs],
where the number of orbits of G is |Orb(G) | = s.
If u and v belong to the same orbit of Gn, then f (Gn, u, x) = f (Gn, v, x). We provide in Table 1 a number of different growth functions in Gn for n = 1,2,3,4. Additionally, we provide the number of vertices and edges, degree set and set of number of vertices with the same degree.
n	|V (Gn) |	|E(G„)|	D(Gn)	VD(Gn)	|Orb(Gn)|
1	4	4	{2}	{4}	1
2	12	16	{2, 4}	{8,4}	2
3	44	64	{2,4, 8}	{32, 8,4}	3
4	172	256	{2,4, 8,16}	{128, 32, 8,4}	5
Table 1: Number of vertices and edges, ascending ordered degree set, set of number of vertices with the same degree and number of orbits in diamond lattice Gn for n = 1, 2, 3,4.
It is not hard to see that
|V (Gn)|
|E(Gn)| D(Gn)
VD(Gn)
2
= 3(2 + 4n)
= 4n
= {2^ = 1, 2,...,n}
(n-i)+1 |i = 1, 2,
, n — 1
} U {22} .
S. J. Sirca and M. Omladic: The JLS model with ARMA/GARCH errors
69
If we are given a root r in Gn this defines a root r in the corresponding Dik (Gn). It is not hard to see that, if dGn (r) = 2® for i e {1, 2,..., n}, then dDi(Gn)(r) = dGn+1 (r) = 2i+1.
We now consider the cases of all possible growth functions of a graph Gn (1 < n < 4) when the root r is a single vertex. First we examine root r of graph Gn which produces an identical growth function as main root vertex rm. If r e V(Gn), r = rm and dGn(r) = dGn (rm),then f (Gn,r,x) = gn(x) and(rm)(Gn) = 4. That means only four vertices in graph Gn (n > 1) produce an identical growth function gn(x) (black and white vertices on Figure 1). Next, for n = 2,3,4 the growth functions of a rooted graph Gn at r, where dc2 (r) = 2, da3 (r) = 4, dG4 (r) = 8 (red vertices on Figure 1) are:
f (G2,r,x) = 1 + 2x + 5x2 + 2x3 + 2x4
f (G3, r, x) = 1 + 4x + 2x2 + 12x3 + 5x4 + 8x5 + 2x6 + 8x7 + 2x8 f (G4, r, x) = 1 + 8x + 4x2 + 8x3 + 2x4 + 24x5 + 12x6 + 24x7 + 5x8 + 16x9 +8x10 + 16x11 + 2x12 + 16x13 + 8x14 + 16x15 + 2x16.
Moreover, using the polynomial p(x) = 1 + 3x + 2x2 + 2x3 and function fn (x) defined in (2.2), we can write the following system of equations that determine the growth functions of a rooted graph Gn at r, where dGn (r) = 2n-1:
f (G3,r,x) = f (G2,r,x2)+22xfo(x)p(x2) f (G4,r,x) = f (G3,r,x2)+23xf1(x)p(x4)
f (Gn, r, x) = f (Gn-1,r,x2) + 2n-1xfn-3(x)p(x2"-2).
For n = 3,4, the growth functions of a rooted graph Gn at r, where dG3 (r) = 2 and da4 (r) = 4 (green vertices on Figure 1) are:
f (G3, r, x) = 1 + 2x + 9x2 +4x3 + 10x4 + 3x5 + 8x6 + 3x7 + 4x8 f (G4, r, x) = 1 +4x + 2x2 + 20x3 + 9x4 + 16x5 + 4x6 + 24x7 + 10x8 + 16x9 +3x10 + 16x11 + 8x12 + 16x13 + 3x14 + 16x15 + 4x16, = f (G3, r, x2)+4xq(x2),
where q(x) = 1 + 5x + 4x2 + 6x3 + 4x4 + 4x5 + 4x6 + 4x7.
Finally, we examine the root r of graph G4 with dG4 (r) = 2, which belongs to one of two orbits. We denote the root in the first orbit by r1 (ur1 e E(G4) and dG4 (u) = 16) and in the second orbit by r2 (ur2 e E(G4) and dG4 (u) = 8). Then the growth functions
f (Gn, r®, x) for n = 4 and i = 1,2 are:
f (G4, r1, x) = 1 + 2x +17x2 +8x3 + 18x4 + 5x5 + 16x6 + 7x7 + 20x8 + 5x9
+16x10 + 6x11 + 16x12 + 6x13 + 16x14 + 5x15 + 8x16 f (G4, r2, x) = 1 + 2x + 9x2 +4x3 + 16x4 + 7x5 + 24x6 + 9x7 + 20x8 + 6x9 +16x10 + 5x11 + 16x12 + 5x13 + 16x14 + 8x15 + 8x16.
2.5 Derivation of the JLS model
The structure of financial markets, given by hierarchical diamond lattice, can be described as follows. Start with a two linked traders. Secondly, replace this link by a diamond, where
70
Ars Math. Contemp. 13 (2017) 125-136
the original traders occupy two diametrically opposite vertices, and the two other vertices are occupied by two new traders. Thirdly, for each of these four links, replace it by a diamond in the same way. If we repeat this process, we get a hierarchical diamond lattice and after n iterations we have § (2 + 4n) traders and 4n links among them.
Johansen et al. [14] constructed so-called JLS model, a non-linear model suitable for the prediction of crash time in both microscopic and macroscopic modelling. This model characterises the occurrence of a crash by the crash hazard rate h (t). All traders are organised into a network and they locally influence each other through this model. The model assumes that agents tend to imitate the opinions of their nearest neighbours. The imitation process is described by the crash hazard rate h (t) with a power law, i.e. dh/dt = Chs, where C is a positive constant and S > 1.
In the JLS model [14] we consider a network of traders, where each trader is indexed by an integer number i = 1,... ,N and N (i) denotes the set of traders who are directly connected to trader i in the network (hierarchical diamond lattice). For simplicity we consider a special case of the Potts model for q = 2. We assume that trader i can be in only one of two possible states at time t: the buy (s¿ = +1) or the sell (s¿ = -1) decision. Then the state of trader i is determined by the following Markov process:
si = sign ( K sj + aei + G ^ jew(i)
where the sign function, sign (x), is equal to +1 (-1) if x > 0 (x < 0), K is positive constant, ei ~ N (0,1) is an independent and identically distributed random variable and term G is a measure of some external influence, which tends to favor state +1 (-1) if G > 0 (G < 0). In this model, K governs the tendency of imitation among traders, while a governs their idiosyncratic behaviour. If we define the average state of the market as M = (1/1) J2i=i si, the susceptibility of the system is defined as x = dE}e ] and measures the sensitivity of the average state of the system to a small external influence.
Further the JLS model assumes that the crash hazard rate behaves in a similar way as the susceptibility in the neighbourhood of a critical point. By considering a hierarchical diamond lattice for the financial market, the dynamics of the crash hazard rate can be described as follows:
h (t) « B0 (tc - t)m-1 + C0 (tc - t)m-1 cos [w ln (tc - t) - ,
where B0, C0 and w are real number.
The dynamics of the price is described as dp = ^ (t) p (t) dt - Kp (t) dj, where p (t) is the price, ^ (t) is time-varying drift and j is a jump process, such that dj = 0 before crash and dj = 1 after the crash occurs at critical time tc. The parameter k determines the loss amplitude associated with the occurrence of a crash.
One assumption of this model is that the price p(t) follows a martingale process, i.e. Et [p(t )] = p (t), Vt > t, where Et [•] represents the conditional expectation on all information available up to time t. Then we have Et [dp] = 0. The dynamics of the jumps is governed by a crash hazard rate h (t) and Et [dj] = h (t) dt. Furthermore, the drift ^ (t) is chosen so that the martingale condition is satisfied, which yields ^ (t) = Kh (t).
Then the simplest form of a log-price dynamics up to the end of the financial bubble can be written as follows:
ln (pt) = A + B (tc - t)m + C (tc - t)m cos (w ln (tc - t) - + ut,
(2.5)
S. J. Sirca and M. Omladic: The JLS model with ARMA/GARCH errors
71
where pt is the price of a stock index or some other specific asset at time t, A = ln (ptc) > 0 is the logarithm of the price at the critical time tc and B < 0 for a growing bubble before the crash. The critical time tc is the end of a bubble and indicates a change to a new market regime, which could be a large crash or a change in the price growth rate. Note that tc is the most probable time for the crash, but there also exists a nonzero probability that the bubble ends without a crash. If C = 0, then the presence of log-periodic behaviour is indicated. The exponent m lies between 0 and 1 to ensure a finite price, even at tc. Parameter w is the frequency of oscillations during the bubble period, while ^ is a phase parameter and lies between 0 and 2n. The error term ut has a zero mean and some standard deviation.
2.6 Estimation (optimisation) of the JLS model
The JLS model in equation (2.5) is described by three linear parameters, A, B and C, and four nonlinear parameters, tc, m, w and For simplicity, we denote the 7-dimensional vector of these parameters by 0 = [A, B, tc, m, C, w, We need to find the vector 0 that best fits the observed log-price time series {ln (p^}^ within the estimation time period [ti, tn], where tn < tc. Although different measures can be used, the most common approach is to minimise the sum of the squared residuals:
where ut is the error term in the JLS model.
The minimisation of such an objective function S (0) is not an easy task due to the presence of many local minima. Therefore, we first reduce the number of parameters by expressing three linear parameters as a function of the remaining four nonlinear parameters, as proposed in [14]. If we rewrite the equation (2.5) as ln (pt) « A + Bft + Cgt, where ft = (tc - t)m and gt = (tc - t)m cos (w ln (tc - t) - then we obtain the estimates of linear parameters A, B and C by using the ordinary least squares (OLS) method:
If we rewrite the system of equations (2.7) in matrix form as XTy = (XTX)^, where
then the well-known solution is / = (XTX)-1XTy.
By reducing the number of parameters from seven to four we decreased the complexity of the optimisation problem. However, we still need to find the global minimum in a 4-dimensional space of the objective function:
(2.6)
(2.7)
Si(0i) = min S(0),
A,B,C
(2.8)
where 01 = [tc, m, w, denotes the vector of nonlinear parameters and S(0) is given in equation (2.6).
72
Ars Math. Contemp. 13 (2017) 125-136
Many different optimisation algorithms have been proposed to estimate the JLS model [30, 14, 11, 27, 8, 18]. In this paper, we use the Differential Evolution (DE) algorithm, which is a simple and efficient population-based search heuristic developed by Storn and Price [33]. We employ the DEoptim function in the R package DEoptim (see [21] for more details) together with the following intervals for each parameter's optimisation:
where t\ and tn are first and last data point in the estimation time period. These intervals are similar to those used in [37].
If the parameter m is too small, then we obtain a bubble with a sudden acceleration at the end, while too large an m corresponds to an almost linear non-accelerating bubble. Similarly, if the frequency w of the oscillations is too small, then the log-periodic oscillations will be too slow, while for too large a value of w, the oscillations are too fast. Therefore, after the estimation of the JLS model, we accept all results that satisfy the following four additional constraints:
B < 0, 0.1 < m < 0.8, 4 < w < 15 and tc < 2tn - tx. (2.11)
Similar bounds for m and w were used in [29]. It is also reasonable to eliminate the results, for which the upper bound of the search space for tc under (2.9) is attained. Therefore we use the last constraint under (2.11). For the sake of brevity, we refer to conditions (2.11) as the LPPL conditions.
2.7 Verification of the JLS model
Important step in building a model is determination of its quality. If the model specification is appropriate, then the residuals should behave like true stochastic components. If this component is white noise, then the residuals should behave like independent (normal) random variables, with a zero mean and some standard deviation.
To investigate the stationarity of the residuals ut in the JLS model (2.5), we perform Phillips-Perron (PP) and Augmented Dickey-Fuller (ADF) unit root tests, where the null hypothesis is the presence of a unit root. Additionally, the Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test is employed for testing the null hypothesis, which posits that the observable time series is stationary. We perform the Shapiro-Wilk test on the residuals as a formal test of normality.
We are also interested in searching for possible dependencies in the JLS model residuals. We apply a runs test to verify the independence of the residuals. Using the autocorrelation function (ACF) of the residuals and squared residuals, we can observe the presence of linear and nonlinear dependence in the residuals. However, for a mixed autoregressive moving average (ARMA) model, it is usually difficult to identify the ARMA orders p and q from these plots. Tsay and Tiao [34] proposed an extended autocorrelation function (EACF) method for model identification that is able to identify mixed ARMA(p, q) models, as well as pure AR(p) and MA(q) models. Details about the EACF method and our proposed algorithm for ARMA order determination can be found in supplementary file (section 1).
We first apply the EACF method on the JLS model residuals to select the appropriate ARMA model. Then we compare this result with two of the information criteria,
tc G [tn + 1, 2tn - ti + 1] , m G [10-5,1 - 10-5] w G [10-5, 40] , ^ G [10-5, 2n - 10-5]
(2.9) (2.10)
S. J. Sirca and M. Omladic: The JLS model with ARMA/GARCH errors
73
namely Akaike's Information Criterion (AIC), given by AIC = log (a2) + 2k, and Schwarz Bayesian Information Criterion (BIC), given by BIC = log (a2 ) + n log (n), where k is the number of estimated parameters, a2 is the estimated error variance of the fitted model, and n is the number of observations. We estimate the parameters of all ARMA(p, q) models with 0 < p,q < 5 and select the ARMA orders by minimising AIC and BIC, respectively. For this purpose, we employ an arima function in the R package TSA.
3 The JLS model with ARMA/GARCH errors
Gazola et al. [8] proposed a log-periodic-AR(1)-GARCH(1,1) model to capture the structure of the JLS model residuals. According to the analysis of residuals described in subsection 2.7, our empirical results on different stock indices show similar results, namely, that the residuals are strongly autocorrelated and in many cases also heteroscedastic. Furthermore, by using the EACF method we also found that an ARMA model is sometimes more appropriate than a pure AR model. Consequently, to capture the behaviour of the error term
ut in (2.5), we propose the following ARMA(p, q)/GARCH(P, Q) error model:
p	q
ut = E Pi Ut-i + nt + E^j nt-j,	(3.1)
i=i j=i
nt = atet,	(3.2)
p	Q
o\ = ao + J2 aknt-k + E ^lat-l,	(33)
k=1 1=1
where et is independent and identically distributed process with a zero mean and unit variance. If q = 0 in equation (3.1), we have a pure AR(p) process and if p = 0 we have a pure MA(q) process. We also include the conditional variance a2, which evolves according to a GARCH(P, Q) process described by equation (3.3). If P = Q =0 in equation (3.3), we have a constant variance and obtain a pure ARMA(p, q) error model.
3.1 Estimation of the JLS model with ARMA/GARCH errors
Let us denote by 0 = [A, B, tc, m, C, w, pi, Qj] the vector of the parameters of the JLS model with ARMA(p, q) errors and by 01 = [A, B, tc, m, C, w, pi, Qj,ak, $], the vector of the parameters of the JLS model with ARMA(p, q)/GARCH(P, Q) errors, where i = 1,..., p, j = 1,..., q, k = 0,..., P, and l = 1,..., Q. In our empirical analysis, we use the following three-step procedure for model identification and estimation of the JLS model with ARMA/GARCH errors:
1.	Estimate the JLS model (2.5), as described in subsection 2.6 and verify that the LPPL conditions under (2.11) are satisfied.
2.	Identify the ARMA orders with EACF method (compare with AIC and BIC) as described in subsection 2.7 and estimate the JLS model with ARMA(p, q) errors. Using the conditional maximum likelihood method, under the normality assumption for et in equation (3.2), the estimates of parameters 0 are obtained by the maximization of the log-likelihood function:
L (0) = -^ log(2n) - ^ log(a2) - £ , (3.4)
t=p+i
74
Ars Math. Contemp. 13 (2017) 125-136
where rjt and a2 = a2 are obtained with equation (3.2), and n is the length of the residuals series.
Perform the Engle's Lagrange Multiplier (LM) test (for lags 10, 12 and 20) to verify the presence of autoregressive conditional heteroscedasticity (i.e. ARCH effect) in the residuals from the estimated model in the previous step. If necessary, identify the GARCH orders with EACF method (compare with AIC and BIC) on the squared and absolute residuals. The maximum likelihood estimation for the JLS model with ARMA(p, q)/GARCH(P, Q) errors can be carried out by maximizing the log-likelihood function:
1 n	1 n 2
L (0i) = -n-p log(2n) - 2 £ log(o?) - 2 £ , (3.5)
t> o a2
--2 - U +
t = p+1 t = p+1 t
where n2 and are obtained with equation (3.2), and n is the length of the residuals series.
For simplicity, if there is an ARCH effect in the residuals, we incorporate a simple GARCH(1,1) model. That means, we are not interested in particular GARCH order, but only in the existence of GARCH structure in the residuals. Therefore, we propose new algorithm for determination if there exist an ARCH effect; see section 1 in supplementary file. Some additional comments on described three-step procedure can be found in supplementary file (section 2).
3.2 Verification of the JLS model with ARMA/GARCH errors
If our proposed model specification is appropriate, than the standardized residuals et = nt/Ut are approximately independent and identically distributed. Therefore, we proceed with the analysis of the standardized residuals from the fitted JLS model with ARMA errors (the model estimated in step 2 of three-step procedure described in previous subsection), if the results in step 3 show that there is no need to incorporate a GARCH process. Otherwise we perform the analysis of the standardized residuals from the fitted JLS model with ARMA/GARCH errors. The standardized residuals are investigated using the same set of tests as described in subsection 2.7.
4 Empirical results 4.1 The data
The dataset consists of daily closing prices for 33 major stock indices worldwide, namely 10 American, 12 European and 11 Asian/Pacific indices. Detailed information about the dataset can be found in supplementary file (section 3). The data were downloaded from Yahoo! finance to the 28th July 2014. According to the available data, we selected stock indices that represent trading activities on main stock exchanges in each geographic region to cover. We did so to cover adequate volume of trading activity in order to assure the generality of our results.
4.2 The JLS model estimation
We apply the procedure as described in section 3 to estimate the JLS model parameters. For all stock indices in the dataset, we generate a set of time windows, each consisting of
S. J. Sirca and M. Omladic: The JLS model with ARMA/GARCH errors
75
500 successive trading days (which is approximately two calendar years). Note, the JLS model is usually estimated on 2-3 years large time window. Each set is obtained by rolling such time windows over the whole dataset with increment of 1 day.
To summarize, the results of the JLS model application to rolling time windows for which the LPPL conditions, as specified in (2.11), are satisfied show, that fraction of such time windows varies between 15.14% and 2.39%. We note that such variation among stock indices reflects also different starting points of the dataset and therefore different bubbles representation across stock indices. To increase the possibility that the selected samples (for further analysis) resemble the bubble periods, we set the minimum number of successive time windows (that satisfy LPPL conditions) to 50. For more details, see section 4 in supplementary file.
4.3 The JLS model residual analysis
In this subsection we investigate the JLS model residuals ut in (2.5) and estimate the JLS model with ARMA/GARCH errors. We perform analysis of 30 stock indices (121 selected samples). For three indices no sample is selected. Detailed information about the selected samples can be found in supplementary file (section 5).
To summarize, results shows that PP, ADF and KPSS tests indicate the stationarity for most cases at the 5% level. For each stock index, we examine any possible dependence in the JLS model residuals by performing a runs test on all time windows of selected samples. For all cases we get p-values of < 0.05, therefore we can reject the null hypothesis of independence at the 5% level. The Shapiro-Wilk test rejects the normality assumption in most cases at the 5% level. For more details, see section 6 in supplementary file.
To see how Gazola et al. [8] proposed model specification holds in the case of our selected data, we first investigate the presence of ARMA orders, different from ARMA(1,0) in selected samples. According to the EACF method the results indicate that the fraction of ARMA(1,0) models varies significantly across the stock indices. We found that, contrary to findings in [8], for only four stock indices the EACF method suggest only the ARMA(1,0) model, while for two indices it suggests only ARMA models with p, q > 0. In most cases the BIC method confirms the results of EACF method, while in general the selected ARMA orders with AIC method are greater compared to those selected by BIC and EACF methods. For more details, see section 7 in supplementary file.
Next, we investigate the existence of GARCH structure in the residuals. We use the following strategy to decide, whether there is no need to incorporate GARCH process: if at least one EACF result on squared or absolute residuals suggest that P = Q = 0 and LM test for at least two lags (10, 12 or 20) confirms that results, or LM test results cannot reject the null hypothesis at the 5% level for all three lags, then the JLS model with only ARMA errors is selected. The results show that, contrary to findings in [8], only for about third of all analysed stock indices the incorporation of GARCH is always necessary and also that the biggest proportion of such indices is in Asian/Pacific region. Here we note that the results comparison between different geographic regions is limited due to different historical data availability. For more details, see section 7 in supplementary file.
As a last step, we proceed with the analysis of the standardized residuals from our proposed generalized JLS (GJLS) model. The Shapiro-Wilk test on standardized residuals of the GJLS model shows quite similar result as JLS model residual analysis. There exists only one index for which the fraction of rejected normality assumption is smaller than 85%. Similarly, in most cases, the PP, ADF and KPSS indicate the stationarity of the
76
Ars Math. Contemp. 13 (2017) 125-136
standardized residuals at the 5% level. The key difference between analysis of the JLS and GJLS model residuals is obtained using a runs test. For 16 indices, performing this test yields fraction smaller than 5% of the rejected null hypothesis at the 5% level. Note that for almost all indices, where this fraction is larger than 5%, using the GJLS model we identify the false JLS model bubble identification. Detailed explanation of the false JLS model bubble identification can be found in supplementary file (section 9).
4.4 The JLS model versus GJLS model
In this subsection we compare the results of the JLS and GJLS model. In doing so, we focus on the parameter of our key interest, critical time tc. For all 121 selected samples, we compute two location parameters for the distribution of tc, namely the mean and the median, and 25-75% quantile interval for tc. For some explanatory comments on the estimation methodology; see section 8 in supplementary file.
Our main results show that GJLS model in general performs better in predicting the actual local or global market peak, denoted by tc obs, which is consistent with findings in [8], that the GJLS model outperforms the JLS model in predicting tc. For the American stock indices the mean or median estimate of the GJLS model is 26 times closer to tc obs, while the JLS model only 14 times. Similar results are also obtained for the European and Asian/Pacific stock indices; see section 9 in supplementary file.
Estimated 25-75% quantile intervals (QI) for tc (calculated with the JLS or GJLS model) cover the tc obs for more than half (sub)samples for stock indices in all three geographic regions. The results differ, however, when comparing the JLS and GJLS models; in the case of the American stock indices the coverage for the JLS model is 14 times, while for the GJLS 18 times. There are also 6 cases for each model, where the estimated QI misses the tc obs for less than 11 trading days. For the European stock indices the QI coverage is 17 (25) times for the JLS (GJLS) model and the number of estimates outside the interval by less or equal to 2 weeks is 11 (12). For the Asian/Pacific stock indices the JLS (GJLS) model QI covers the tc obs 5 (10) times, with 3 estimates outside the intervals for less or equal to 2 weeks for each model. For more details, see section 9 in supplementary file.
Comparing the mean and median estimates across the JLS and GJLS models for the American stock indices, the former performs better in 24 cases, while the later in 14 cases. Similar outcome is observed also for the European stock indices. For Asian/Pacific indices the mean and median perform almost equally good, but we note that for this region we have considerably lower number of (sub)samples. It therefore makes sense to take into account also the choice of the point estimate measure. For more details, see section 9 in supplementary file.
5 Conclusions
In this paper we consider the Johansen-Ledoit-Sornette (JLS) model, which describes the behaviour of stock index prices during a bubble regime and also predicts the most probable time for a change in the regime. We introduce the diamond lattice operation Di which is a composite operation of parallelization and one-dimensional subdivision, to describe the construction of the hierarchical diamond lattice used in JLS model derivation. It turns out that the operation Di has interesting properties that we are investigating. The results will be published elsewhere.
The idea of our JLS model generalization is motivated by the behaviour of the JLS
S. J. Sirca and M. Omladic: The JLS model with ARMA/GARCH errors
77
model residuals. The results of our analysis reveal the presence of a strong autocorrelation in JLS model residuals and heteroscedasticity in many cases. To incorporate these two properties into the model, we propose an ARMA(p, q)/GARCH(P, Q) error model and a methodology for model identification, parameter estimation and its verification.
As the first part of our analysis we investigate the behaviour of the JLS model residuals. To assure the generality of the results, we perform an analysis over a large size of samples for 33 stock indices from three broader world geographic regions. Our results suggest that there is no general rule for determination of ARMA/GARCH specification of the JLS model as the proportion of ARMA models, excluding ARMA(1,0) and proportions of ARMA models without GARCH specification varies significantly across the samples and also within each geographic region. We take this result as a reminder for careful ARMA/GARCH order determination, when analysing particular time period for stock indices.
In the second part of the analysis we show that specified JLS model generalization outperforms the JLS model estimates of critical time tc. Using the mean value and the median of tc parameter estimates, the results show smaller deviations from the actual crash dates for the GJLS model. Moreover, we also show that 25-75% quantile intervals more often cover the tc parameter estimates of the GJLS model than for the JLS model. Further more, with GJLS model we are able to detect false alarms, meaning that there is no actual bubble period, which is otherwise identified with JLS model.
Acknowledgements
This research was funded in part by the European Union, European Social Fund, Operational Programme for Human Resources Development for the Period 2007-2013 under grant P-MR-08/63; and Slovenian Research Agency under grant P1-0222.
References
[1]	A. N. Berker and S. Ostlund, Renormalisation-group calculations of finite systems: order parameter and specific heat for epitaxial ordering, J. Phys. C 12 (1979), 4961.
[2]	H.-C. G. v. Bothmer et al., Significance of log-periodic signatures in cumulative noise, Quant. Finance 3 (2003), 370-375.
[3]	A. Clark, Evidence of log-periodicity in corporate bond spreads, Phys. A 338 (2004), 585-595.
[4]	M. del Rio Francos, Chamfering operation on k-orbit maps, Ars Math. Contemp. 7 (2014), 507-524.
[5]	B. Derrida, L. De Seze and C. Itzykson, Fractal structure of zeros in hierarchical models, J. Statist. Phys. 33 (1983), 559-569, doi:10.1007/BF01018834, http://dx.doi.org/10. 10 0 7/BF01018834.
[6]	J. A. Feigenbaum, A statistical analysis of log-periodic precursors to financial crashes, Quant. Finance 1 (2001), 346-360, doi:10.1088/1469-7688/1/3/306, http://dx.doi.org/10.
1088/1469-7688/1/3/30 6.
[7]	P. Fowler and T. Pisanski, Leapfrog transformations and polyhedra of Clar type, J. Chem. Soc., Faraday Trans. 90 (1994), 2865-2871, doi:10.1039/FT9949002865.
[8]	L. Gazola, C. Fernandes, A. Pizzinga and R. Riera, The log-periodic-AR(1)-GARCH(1,1) model for financial crashes, Eur. Phys. J. B 61 (2008), 355-362.
78
Ars Math. Contemp. 13 (2017) 125-136
[9] P. Gnaciñski and D. Makowiec, Another type of log-periodic oscillations on Polish stock market, Phys. A 344 (2004), 322-325.
[10]	I. Hubard, M. del Rio Francos, A. Orbanic and T. Pisanski, Medial symmetry type graphs, Electron. J. Combin. 20 (2013), Paper 29, 28.
[11]	E. Jacobsson, How to predict crashes in financial markets with the Log-Periodic Power Law, Master's thesis, Department of Mathematical Statistics, Stockholm University, 2009.
[12]	Z.-Q. Jiang, W.-X. Zhou, D. Sornette, R. Woodard, K. Bastiaensen and P. Cauwels, Bubble diagnosis and prediction of the 2005-2007 and 2008-2009 Chinese stock market bubbles, J. Econ. Behav. Organ. 74 (2010), 149-162.
[13]	A. Johansen, Discrete scale invariance and other cooperative phenomena in spatially extended systems with threshold dynamics, Ph.D. thesis, Niels Bohr Institute, University of Copenhagen, 1997.
[14]	A. Johansen, O. Ledoit and D. Sornette, Crashes as critical points, Int. J. Theor. Appl. Finance 3 (2000), 219-255.
[15]	A. Johansen and D. Sornette, Financial "anti-bubbles": Log-periodicity in gold and Nikkei collapses, Int. J. Mod. Phys. C 10 (1999), 563-575.
[16]	A. Johansen and D. Sornette, Log-periodic power law bubbles in Latin-American and Asian markets and correlated anti-bubbles in Western stock markets: An empirical study, Int. J. Theor. Appl. Finance 4 (2001), 853-920.
[17]	M. D. Konig and S. Battiston, From graph theory to models of economic networks. A tutorial, in: Networks, Topology and Dynamics, Springer, 2009, doi:10.1007/978-3-540-68409-1_2, http://dx.doi.org/10.100 7/97 8-3-540-6840 9-1_2.
[18]	V. Liberatore, Computational LPPL fit to financial bubbles, arXiv (2011), accessed 30 January 2 011.
[19]	L. Lin, R. Ren and D. Sornette, A consistent model of 'explosive' financial bubbles with mean-reversing residuals, Swiss Finance Institute Research Paper (2009).
[20]	R. Matsushita, S. Da Silva, A. Figueiredo and I. Gleria, Log-periodic crashes revisited, Phys. A 364 (2006), 331-335.
[21]	K. Mullen, D. Ardia, D. L. Gil, D. Windover and J. Cline, DEoptim: An R package for global optimization by differential evolution, J. Stat. Softw. 40 (2011), 1-26.
[22]	A. Orbanic, D. Pellicer, T. Pisanski and T. W. Tucker, Edge-transitive maps of low genus, Ars Math. Contemp. 4 (2011), 385-402.
[23]	A. Orbanic, D. Pellicer and A. I. Weiss, Map operations and k-orbit maps, J. Combin. Theory Ser. A 117 (2010), 411-429, doi:10.1016/j.jcta.2009.09.001, http://dx.doi.org/10. 1016/j.jcta.2009.09.001.
[24]	T. Pisanski and M. Randic, Bridges between geometry and graph theory, in: Geometry at work, Math. Assoc. America, Washington, DC, volume 53 of MAA Notes, pp. 174-194, 2000.
[25]	T. Pisanski and T. W. Tucker, Growth in repeated truncations of maps, Atti Sem. Mat. Fis. Univ. Modena 49 (2001), 167-176, dedicated to the memory of Professor M. Pezzana (Italian).
[26]	T. Pisanski and A. Zitnik, Representing graphs and maps, in: Topics in topological graph theory, Cambridge Univ. Press, Cambridge, volume 128 of Encyclopedia Math. Appl., pp. 151-180, 2009.
[27]	X. Qin and T. G. K. Randolph, Bubbles, can we spot them? Crashes, can we predict them?, PhD, Vienna Graduate School of Finance (2005), http://repec.org/sce2 0 0 5/up. 18847.1107074423.
S. J. Sirca and M. Omladic: The JLS model with ARMA/GARCH errors
79
[28]	A. D. Sokal, Chromatic polynomials, Potts models and all that, Phys. A 279 (2000), 324-332.
[29]	D. Sornette and A. Johansen, Significance of log-periodic precursors to financial crashes, Quant. Finance 1 (2001), 452-471.
[30]	D. Sornette, A. Johansen and J.-P. Bouchaud, Stock market crashes, precursors and replicas, J. Phys. I 6 (1996), 167-175.
[31]	D. Sornette, R. Woodard and W.-X. Zhou, The 2006-2008 oil bubble: Evidence of speculation and prediction, Phys. A 388 (2009), 1571-1576.
[32]	D. Sornette and W.-X. Zhou, The US 2000-2002 market descent: How much longer and deeper?, Quant. Finance 2 (2002), 468-481.
[33]	R. Storn and K. Price, Differential evolution—a simple and efficient heuristic for global optimization over continuous spaces, J. Global Optim. 11 (1997), 341-359, doi:10.1023/A: 1008202821328, http://dx.doi.org/10.1023/A:1008202821328.
[34]	R. S. Tsay and G. C. Tiao, Consistent estimates of autoregressive parameters and extended sample autocorrelation function for stationary and nonstationary ARMA models, J. Amer. Statist. Assoc. 79 (1984), 84-96, http://links.jstor.org/sici?sici= 0162-1459(198403)79:385<84:CEOAPA>2.0.CO;2-V&origin=MSN.
[35]	D. J. A. Welsh and C. Merino, The Potts model and the Tutte polynomial, J. Math. Phys. 41 (2000), 1127-1152, doi:10.1063/1.533181, probabilistic techniques in equilibrium and nonequilibrium statistical physics, http://dx.doi.org/10.1063Z1.533181.
[36]	J. H. Wosnitza and C. Denz, Liquidity crisis detection: An application of log-periodic power law structures to default prediction, Phys. A 392 (2013), 3666-3681, doi:10.1016/j.physa.2013. 04.009, http://dx.doi.Org/10.1016/j.physa.2013.04.009.
[37]	W. Yan, R. Woodard and D. Sornette, Diagnosis and prediction of market rebounds in financial markets, Swiss Finance Institute Research Paper (2010).
[38]	W. Yan, R. Woodard and D. Sornette, Leverage bubble, Phys. A 391 (2012), 180-186.
[39]	W.-X. Zhou and D. Sornette, Statistical significance of periodicity and log-periodicity with heavy-tailed correlated noise, Int. J. Mod. Phys. C 13 (2002), 137-169.
[40]	W.-X. Zhou and D. Sornette, Evidence of a worldwide stock market log-periodic anti-bubble since mid-2000, Phys. A 330 (2003), 543-583, doi:10.1016/j.physa.2002.12.001, http:// dx.doi.org/10.1016/j.physa.2002.12.001.
[41]	W.-X. Zhou and D. Sornette, Nonparametric analyses of log-periodic precursors to financial crashes, Int. J. Mod. Phys. C14 (2003), 1107-1125, doi:10.1142/S0129183103005212, http: //www.worldscientific.com/doi/abs/10.1142/S0129183103005212.
[42]	W.-X. Zhou and D. Sornette, Is there a real-estate bubble in the US?, Phys. A 361 (2006), 297-308.
[43]	W.-X. Zhou and D. Sornette, Analysis of the real estate market in Las Vegas: Bubble, seasonal patterns, and prediction of the CSW indices, Phys. A 387 (2008), 243 - 260, doi:http://dx. doi.org/10.1016/j.physa.2007.08.059, http://www.sciencedirect.com/science/ article/pii/S0378437107008473.
[44]	W.-X. Zhou and D. Sornette, A case study of speculative financial bubbles in the South African stock market 2003—2006, Phys. A 388 (2009), 869 - 880, doi:http://dx.doi.org/10.1016/j. physa.2008.11.041, http://www.sciencedirect.com/science/article/pii/ S0378437108009709.