Fractional virus epidemic model on financial networks

Abstract In this study, we present an epidemic model that characterizes the behavior of a financial network of globally operating stock markets. Since the long time series have a global memory effect, we represent our model by using the fractional calculus. This model operates on a network, where vertices are the stock markets and edges are constructed by the correlation distances. Thereafter, we find an analytical solution to commensurate system and use the well-known differential transform method to obtain the solution of incommensurate system of fractional differential equations. Our findings are confirmed and complemented by the data set of the relevant stock markets between 2006 and 2016. Rather than the hypothetical values, we use the Hurst Exponent of each time series to approximate the fraction size and graph theoretical concepts to obtain the variables.

variable population size and a delay [14] have been studied by several researchers. In several studies on computer virus epidemiology [15][16][17], continuous state models are used to represent degradation. This approach has numerous advantages to two-state models such as SIS [18].
In this study, we consider the economic crisis or stress as a virus contagious on a financial network of globally operating stock markets, and represent it by the fractional order differential equation system as Since fractional order derivatives emerge as a powerful tool to describe the memory and hereditary properties, they are frequently used in various areas such as fluid mechanics, biology, physics, and economics. Fractional calculus is an effective way of incorporating memory effects. The kernel of power-law defining the fractional relaxation equation presents a long-term memory [19,20]. For the analytical and numeric solutions of fractional order differential equations and their systems, numerous methods such as Laplace transform method [21][22][23], iteration methods [24,25], differential transform method [26,27], Adomian decomposition method [28], and the Fourier transform method [29,30] are proposed by several researchers. The paper is organised as follows: In Section 2, we give the basics of the graph theory to construct stock market network and preliminaries about the fractional calculus. For more details on the graph theory we refer readers to [31]. In Section 3, we first construct the network with respect to correlation distance of each stock market's daily logarithmic return. For the edge forming, we introduce a method to determine threshold value. The model (1) can be studied in two ways: commensurate and incommensurate systems. In the case of the equality of the fraction sizes we obtain the commensurate system, and the case of the inequality of different fraction sizes yields the incommensurate system. For the commensurate system we use an analytical method that is based on Mittag-Leffler functions applied for the fraction size˛D 0:5 in Section 4. By choosing˛D 0:5, we discard the autocorrelation of each time series. Again in Section 4, we apply the differential transform method for the solution of the incommensurate system. For the commensurate system, we determine the fraction sizes with respect to Hurst exponent of each time series. The concluding remarks on this study are given in Section 5.

Preliminaries
In this section, we present mathematical preliminaries and concepts of fractional calculus and graph theory that we used throughout the paper. For constructing the financial network we use undirected simple graphs and some of their classes, and spectrum to determine threshold boundaries. For the epidemic model acting on a network, we use the Caputo fractional derivatives. Throughout the paper represents the Gamma function which is an extension of the factorial function.

Fractional calculus
Fractional calculus generalizes the notions of integer-order differentiation and integration to the fractional order [22,32]. Here, we give some basic definitions and properties of fractional calculus theory that will be used throughout the paper.
Definition 2.1 ( [33]). For f .x/ 2 C.a; b/ and n 1 <˛Ä n, we let There are several generalization of the fractional derivatives. The most common ones are the Riemann-Liouville and Caputo senses. When solving fractional differential equations, the Caputo's definition does not require fractional order initial conditions [22]. We have chosen to use the Caputo fractional derivative because it allows traditional initial conditions to be included the formulation of the problem. The fractional differentiation in the Riemann-Liouville sense is defined by for n 1 Ä˛< n. The relation between the Riemann-Liouville and Caputo operators can be given by for n 2 Z C [22]. For the sake of simplicity, we denote the Caputo operator as C a Dt D Dą throughout the study. Since our formulation involves only the initial conditions we choose a D 0.
As the solutions of integer order differential equations are expressed in the terms of exponential functions, the solutions of fractional order differential equations can be expressed in the terms of Mittag-Leffler functions. More details on Mittag-Leffler functions can be found in [22,34] and the references therein.

Graph theory
One of the most powerful tools to represent real world problems which involve relationships is a simple graph. A simple graph is the tuple G D .V; E/ where V is the set of vertices (or the nodes representing the agents in interaction) and E is the set of edges (or the links representing the interactions). Elements of E are the unordered pairs of vertices; i.e., E D fe D .v i ; v j /jv i ; v j 2 V g.
A path between any vertices v i and v j can be defined as the sequence of edges whose end points are v i and v j . The graph G D .V; E/ is called connected if there is a path between any vertices. The complete graph is an undirected graph with every pair of apart vertices connected by an edge. A k-clique is a subset of k-vertices of an undirected graph such that its induced subgraph is complete. In the case of real world data representation, each edge in E may be assigned by a non-negative numerical value. This value is called as a weight. The triple .V; E; w/ is called as a weighted graph for the mapping w W E ! R C .
To represent simple and undirected graph, we use a binary matrix whose elements are This matrix is called adjacency matrix and is symmetric by the definition. For a vertex v in an undirected graph G, the number of edges incident to that vertex is called a degree and let us denote it by d v . For the graph G D .V; E/, the Laplacian Matrix of G is the matrix whose entries are The graph Laplacian does not depend on an ordering of the vertices of G. For the spectrum of the Laplacian S G D f 1 ; : : : ; m g, it is well known that all of the eigenvalues are positive semi-definite; i.e. i 0 with the least one 0. For a graph with nonnegative weights, the multiplicity k of the eigenvalue 0 of L G equals the number of connected components in the graph [35].
One of the best known graph problems is to find a tree structure that has minimum weight. There are several approaches to solve this problem and the resulting tree is known as Minimum Spanning Tree (MST) [36][37][38]. MST is also known to bring hierarchical structures in complex networks [39,40]. The latter studies show that the tree like structures that involves cliques serve more efficiently to bring hierarchies [41,42]. Planar graphs have the same hierarchical structure as MST but they contain a larger amount of edges, loops and cliques. The idea of the construction of planar graphs is based on connecting the most correlated agents iteratively while constraining the resulting network to be embedded on a surface with genus g. In [42], authors briefly studied the special case for g D 0; i.e., the graph embedded on a sphere, and called it as the Planar Maximal Filter Graph (PMFG). PFMG can be seen as the topological triangulation of the sphere, henceforth it is only allowed to have three or four cliques [43].

The model
In this study, a financial network of globally operating stock markets are modelled by a simple undirected graph For the daily closure price C l i of the i -th stock market, the daily logarithmic return R i is calculated as and to obtain links between the stock markets, we use the Pearson correlation of each stock market as where < :: > is a temporal average performed on the trading days.
To avoid the non-positive weights on edges, we introduce the distance function involving correlations as C orrDist WD p 2.1 ij /=2. While correlation coefficients vary between 1 to 1, the values of C orrDi st vary 0 to 1.
To determine the links, we use the following formation rule: Here, T hV is an empirical threshold value. Our method to determine the network initially starts with a complete graph. The spectrum of this complete graph involves only one 0 eigenvalue. Then, we subdivide OE0; 1 closed interval with 1= h step. The natural candidate to threshold value is i= h for i D 0; : : : ; h. Hence the lower boundary for threshold value is 0; i.e., the initial complete graph. At certain point between 0 and 1, the graphs formed by the forementioned formation rule will have more than one component. To check this, counting the 0 eigenvalues of each graph spectrum will reduce the computational complexity.
In this study, we let h D 100 and determined the boundaries for T hV as OE0; 0:86. The multiplicity of the 0 eigenvalues of each graph for i=100; i D 0; : : : ; 100 is given in Figure 1.
To determine the convenient T hV in boundary, we need to introduce a concept called the disparity measure. For each graph with respect to T hV 2 OE0; 0:86, it is possible to obtain PMFGs. An analysis on all of the 4-cliques in the PMFG reveals a high degree of homogeneity with respect to the 17 globally operating stock markets. The mean of disparity measure < y > can be defined analogously to [42] as the mean of over the clique, where i is a generic element of the clique and C orrDi st .i; j /: It should be noted that this measure is only valid for the connected graphs [42], henceforth we computed it only for T vH 2 OE0; 86. More interestingly, local boundaries occur for the mean disparity measure. For T hV 2 OE0:81; 0:86, the formed graphs yield infinite disparity measures since the s i values tend to 0. Therefore, we choose the T hV as 0:8 to model the network more efficiently. In Figure 2, the values of the disparity measure for the T hV 2 OE0; 0:8 are given. The network under consideration is given in Figure 3 as the triple G D .V; E; w/. Each vertices represent the corresponding stock market. The set of edges are formed with respect to C orrDi st Ä 0:8. The natural weight defined on network is the C orrDi st value. The central stock markets, that have the minimal greatest graph distance to other vertices, are DOW, KOSPI, NIKKEI, and SMI. The rest serve as peripheral nodes. Also DOW, EUROSTOXX, KOSPI and ATX, DOW, MERVAL form 3-cliques, respectively. In our study, the vertices of the network are considered as the infectious components of the virus epidemic model. The spreading of the virus is through the noninfectious components, that are the edges of the network. Hence, the components of the network are assumed to be nonhomogeneous. This infection can be seen as the global economic crisis or stress analogously in financial networks. Once an economic crisis occurs in a stock market, the crisis replicate itself into other stock markets that are the adjacent to infectious one in the network. Progressively, more and more stock markets will be influenced by the crisis exponentially, and then the vertices of whole connected network will be under crisis immediately. In most cases, all vertices would be affected before the recovering process of an individual vertex. Hence, a two-state model is not sufficient to model the crisis phenomena in the financial network.
An advantageous alternative to a two-state model is the continuous-state model where the state of each vertex takes real values to represent the crisis state. The continuous state model is defined on the state space 2 OE0; 1. The "0" state is representing the perfect functioning that is there is no crisis, and "1" is the complete crisis state. i .t/ denotes the state index of vertex i at time t . The fractional process plays key role to capture long-range memory and self-similarity behaviour in complex systems. Therefore, we describe our model for financial crisis spreading as an epidemic fractional differential equation. By the aforementioned assumptions, we set i .t / D 0 if the vertex v i is at healthy state. The fractional deviation from the normal state represents the level of crisis occuring in the vertex. The weight w ij is the strength of the interaction between the vertices v i and v j . As we represent our network with the weighted graph G D .V; E; w/, the weights are computed by w ij D C orrDi st .v i ; v j /. For the vertex v i 2 G, its state is affected by the total impact of neighboring vertices v j 2 G. Hence the ability of the vertex v i to overcome a crisis is proportional to its vertex degree d v i . Taking all these assumptions into consideration, the status change of the vertex v i can be expressed as the following epidemic fractional differential equation: where N i is the set of the neighboring vertices of the vertex v i and the˛is the fractional dimension. In our model we use the Hurst exponent (H ) to determine the fractional dimension empirically. H is a statistical measure of long-term memory of time series. Basically, it relates the autocorrelations of the time series. The measurement can be done by where EOEx is the expected value, R.n/ is the range of the first n values, S.n/ is their standard deviation, C is a constant. The value of H varies on .0; 1/ by the definition [44]. H can be referred to as the index of dependence and quantifies the relative tendency of a time series either to regress strongly to the mean or to a cluster in a direction. If 0:5 < H < 1 the time series has a long-term positive autocorrelation while for 0 < H < 0:5 it has a long-term negative autocorrelation. For H D 0:5 can indicate a completely uncorrelated series [44]. The computed H values of the stock market time series used in this study are given in Table 1.
In this study, we examine two cases for the fractional dimension˛. First case is for˛D 0:5, that is the model discarding the autocorrelation of each stock market. Then, the commensurate linear fractional order system can be expressed as where t 2 , -D OE 1 ; : : : ; N , W is the correlation distance matrix, D v is the diagonal matrix with T r.D v / D P N iD1 d v i , and R 0 R N is the vector whose i -th entry is R i .0/, and N is the number of vertices in the network. The stability of such system is studied in [45] where necessary and sufficient conditions are derived. The following theorem is an analogues to similar results.
Theorem 3.1. The autonomous fractional order linear system (4) is asymptotically stable iff j arg.spec.W D v /j > =2. In this case the components of the state decay towards 0 like t F or the case˛D 1; if no poles of the system (4) lie in the closed right-half plane, then the stability occurs. Hence, this is consistent with the results for ordinary systems [46].
The second case is˛i D H i , that is the fractional rate of change is equal to Hurst exponent of each stock market. In this case our model includes the autocorrelation, positive or negative, into consideration. Then, the incommensurate linear fractional order system can be expressed as In this system, only the matrix D˛differs from the system (4). Here D˛D D˛1 t 0 ; : : : ; D˛N t 0 and D˛i t 0 indicates the Caputo fractional derivative of i of order˛i . For the positive rational numbers˛i , i D 1; : : : ; N , the following stability result is proposed in [47], and similar results can be derived as follows:

Solution
In this section, we give two solutions for the systems (4) and (5). The existence and uniqueness of the solution of such commensurate system (4) is studied in [21] and the analytical solution is presented in the terms of Mittag-Leffler functions. For the incommensurate system (5), the approximate solution is obtained by using a differential transform method. This method is presented in [27] and the solutions are obtained accurately. The matrix W with each entry is the correlation distance of each adjacent stock market and is represented in Figure 4.

Analytical solution to commensurate system
To obtain general solution of the commensurate fractional order linear system (4), we seek solutions of the form where the constant and the vector u need to be determined and E˛is the Mittag-Leffler function. Substituting the equation (6) into the system (4) yields u E˛. t˛/ D .W D v /uE˛. t˛/ and upon cancelling the matrix equation where the I is the N N identity matrix. Hence, the vector u given by equation (6) is the solution to the system (4). For this solution it is provided that is an eigenvalue and u is an associated eigenvector of the matrix D v W. Henceforth, the general solution to commensurate system (4) is given by where c 1 ; : : : ; c n are arbitrary constants, 1 ; : : : ; n are the eigenvalues, and u .1/ ; : : : ; u .n/ are the corresponding eigenvectors.
To determine the coefficients c 1 ; : : : ; c n we use the initial condition matrix By the determination of the unknown terms, the analytical solutions of the system (4) operating on the stock market network presented in Section 3 are given in Figure 5-7. In these figures, horizontal axis represents the time variable t 2 OE0; 1 and vertical axis represents the state function .t /. We need to note that our solution in the form of Equation (9) of the system (4) does not include complex or multiple eigenvalues.
In Figure 5, we first study the state function for t 2 OE0; 0:1. In the very first steps the crisis affects all stock markets. However, some of them resist to crisis at certain point and start to recover their state. Hence we group stock markets as descending ( Figure 5a) and ascending (Figure 5b) states. In Figure 5a; EUROSTOXX, SMI, and TELAVIV are the most vulnerable stock markets in the time of crisis and DOW, NASDAQ, and SP500 are the first stock markets resisting the crisis as seen in In Figure 5b. In Figure 5b, it can also be seen that after KOSPI resist to crisis, it is again affected in OE0:09; 0:1 and loose its healthy state. Another Asian stock market NIKKEI is acting in opposite way; it first get affected by the crisis, then strongly resists. In Figure 6, the similar resistances and effects can be seen. In Figure 6a, AEX and BIST resist to crisis and recover their state. In Figure 6b, DAX, IPC, and TSEC get affected by the crisis while others remain in healthy state. While t tends to 1, the stock markets exponentially get affected or resist crisis as can be seen in Figure 7.

Differential transform solution to incommensurate system
In this section, we apply the differential transform method introduced in [27] to obtain the approximate solutions for the incommensurate system (5). It shall be noted that the fractional degrees are different, hence an analytical solution to the system (5) cannot be obtained. The differential transform method a is numerical method based on the series expansion. If we expand the analytical and continuous function f .t / in terms of a fractional power series, then we obtain here˛is the fractional degree and M.i / is the fractional differential transform of .t / [26].
For the Caputo sense fractional differential equations initial conditions are not restricted to fractional orders. Therefore, the transformation of the initial conditions are defined in [26] as follows: where˛is the order of the fractional differential equation. Let N i be the differential transform of i for i D 1; 2; : : :. The following theorems serve as the fundamentals to the differential transform method we use in this study: where H i is the Hurst exponent of i , w ij is the .i; j /-th entry of the weight matrix, and d i is the i-th diagonal entry of the degree matrix. By considering Theorems 4.1-4.2, the system (11) can be transformed to The solution functions of the systems (12) with the initial conditions (13) can be obtained as follow: The results are shown in Figure 8. In the very first steps the crisis affects all stock markets. With the change in the healthy state, the resistance occurs in two manners, that is high and low. In Figure 8a we show the stock markets with low resistance. These are the stock markets that diverge from their equilibrium states. Amongst them, TSEC has the lowest and KOSPI has the highest low resistance. In Figure 5b we show the stock markets with high resistance. DOW is the only stock market that can rapidly recover its unhealthy state. Besides, SP500 and NASDAQ also resist to crisis strongly. All stock markets represented in Figure 8b tend to seek an equilibrium state.

Conclusions
Usually the real world problems are described by using ordinary differential equations. However, it is a hard task to be completed to tackle these problems whenever the memory effect is included. The memory effect is an efficient concept for the models of complex systems. In this study, we propose a nonlinear fractional order model in order to explain and understand the epidemic behaviour of a global economic crisis. Before introducing the model, we first construct the financial network of globally operating 17 stock markets as an undirected graph. The links are formed by a distance function involving the correlation of each time series of the stock markets operating from 02.01.2006 to 15.02.2016. To determine the boundaries of the threshold, we introduce a method based on the spectrum of the graph. For the threshold value to form links we introduce an analogous disparity measure involving the correlation distances in the cliques. The threshold value for the network is determined as 0.8. Larger partitions would yield more accurate threshold value, however it would also exponentially grew the time complexity. Hence, we determine the partition as h D 100. This model can be studied in two ways. First is the commensurate one that is˛1 D : : : D˛1 7 D˛. The commensurate system is the one discarding the autocorrelation of each time series when˛D 0:5. We need to note that the commensurate system does not include multiple or complex eigenvalues in our network. The second system is the one with different fractional order, that is the incommensurate system. Rather than the hypothetical fractional orders, we use the Hurst exponent which is the measure of the autocorrelation of a time series. To solve incommensurate system, we use the differential transform method. This method is convenient for our study because its time complexity is low and it does not have any restriction on the time domain. By using differential transform method, it is possible to forecast the future states of the model as in Figure 9. Neo-classical economic theory states that markets always attempt to reach their equilibrium state. In the time of crisis, this attempt is made through the adjustments done by the institutions. The commensurate system we present in this study is discarding the autocorrelations. Hence, the solution let us to determine which markets can resist or not in the time of crisis. The details on this resistance is explained in Section 4.1. Since the incommensurate system we present in this study considers the autocorrelations, we assume there is no adjustment by the institutions. In this case, some of the markets can converge to their equilibrium state while the rest diverge. By Section 4.2 it can be concluded that the converging markets are with the negative autocorrelation and are globally strong markets, such as DOW, SP500, NASDAQ, FTSE, and NIKKEI.
It is concluded that the fractional order epidemic models we present in this study produces results that agree very well with the real data. Moreover, these models can provide useful information for the understanding, prediction, and control of the global economic crisis.