# Markov Chains application to the financial-economic time series prediction

Vladimir Soloviev<EMAIL_ADDRESS>Cherkasy National University named
after B. Khmelnitsky, Cherkassy, Ukraine  Vladimir Saptsin
<EMAIL_ADDRESS>Kremenchug National University named after M.
Ostrogradskii, Kremenchuk, Ukraine  Dmitry Chabanenko<EMAIL_ADDRESS>Cherkasy National University named after B. Khmelnitsky, Ukraine

###### Abstract

In this research the technology of complex Markov chains is applied to predict
financial time series. The main distinction of complex or high-order Markov
Chains and simple first-order ones is the existing of aftereffect or memory.
The technology proposes prediction with the hierarchy of time discretization
intervals and splicing procedure for the prediction results at the different
frequency levels to the single prediction output time series. The hierarchy of
time discretizations gives a possibility to use fractal properties of the
given time series to make prediction on the different frequencies of the
series. The prediction results for world’s stock market indices is presented.

Prediction, time series, complex Markov chains, discrete time, fractal
properties, discrete Fourier prediction.

###### Contents

  1. 1 Introduction
  2. 2 Analysis of prominent publications relevant to the subject
  3. 3 Aims of the paper, problem statement
  4. 4 Classical modeling problems of ESE systems dynamics
  5. 5 Modern concepts in ESE systems modeling
  6. 6 Markov chains prediction technology
  7. 7 Prediction construction algorithm
  8. 8 States in complex Markov chains and approaches for defining them
  9. 9 Step-by-step prediction procedure. Defining the most probable state on the next step, prediction scenarios
  10. 10 Time increments hierarchy and splicing procedure
  11. 11 Results of stock indices prediction
  12. 12 Conclusions and further work

##  1 Introduction

Successful modeling and prediction of processes peculiar to complex systems,
such as ecological, social, and economical (ESE) ones, remain one of the most
relevant problems as applied to the whole complex of natural, human and social
sciences (r001SamarskyMihaylov01 ; IvakhnenkoMGUA68 ;
r002BogoboyashyKurbanov04 ; AndersenGluzmanSornette2000 ;
PincakCurrencyStringTheory2011 ; PincakStringPrediction2011 ).

The diversity of possible approaches to modeling such systems and, usually,
more than modest success in the dynamics prediction, compel us to look for the
reasons of failure, finding them not only in details, but also in the
axiomatics, which relates to problem statement, chosen modeling methods,
results interpretation, connections with other scientific directions.

With the appearance of quantum mechanics and relativity theory in early
twentieth century new philosophical ideas on physical values, measuring
procedures and system state have been established, the ones that are
completely different from Newtonian notions r003ElutinKrivchenko76 ;
r004LandauQuantNerelyativistic .

For more than 70 years basic concepts of classical and neoclassical economic
theories have been discussed by leading scientists, generating new approaches
r005SapirEconTheoryNeodnorSyst . The general systems theory has acquired
recognition in the middle of the 20th century giving way to development of the
new, systemic, emergent, and quantum in essence approach to investigation of
complex objects, which postulates the limited nature of any kind of modeling
and is based upon fixed and closed system of axioms r006BertalanfyGenSyst62 .

However, the development of this new philosophical basis of ESE systems
modeling is still accompanied with numerous difficulties, and new principles
are often merely declared.

Current research is devoted to investigation and application of the new
modeling and prediction technology, suggested in r007KurbanovSaptsin07 ;
r008SaptsinMarkovPaper09_ , based on concepts of determined chaos, complex
Markov chains and hierarchic (in terms of time scale) organization of
calculating procedures.

##  2 Analysis of prominent publications relevant to the subject

Prediction of financial-economic time series is an extremely urgent task.
Modern approaches to the problem can be characterized by the following
directions: 1) approximation of a time series using an analytical function and
extrapolation of the derived function towards future – so-called trend models
r009LukashinAdapt03 ; 2) investigation of the possible influence various
factors might have on the index, which is being predicted, as well as
development of econometric or more complicated models using the Group Method
of Data Handling (GMDH) IvakhnenkoMGUA68 ; r010ZaychenkoMonogr08 ; 3) modeling
future prices as the decisions-making results using neuronal networks, genetic
algorithms, fuzzy sets r010ZaychenkoMonogr08 ; r011EzhovShumsky98 ;
r012Zaencev .

Unfortunately, these techniques don’t produce stable forecasts, what can be
explained by complexity of the investigated systems, constant changes in their
structure. Although we are trying to join these directions in one algorithm,
it is the latter option that we prefer, with it consisting in creating a model
adequate to the process generating a price time series r013ChabM10_ . This
very approach gives a chance to approach the complexity of the system, which
generates the observed series, develop the model and use its properties as the
prognosis.

##  3 Aims of the paper, problem statement

Assume the time series is set by a sequence of discrete levels with constant
step of time sampling $\Delta t$. We need to generate variants of the time
series continuation (prognosis scenarios) according to the relations between
the sequences of absolute and relative changes discovered with the help of
complex Markov chains.

##  4 Classical modeling problems of ESE systems dynamics

Another peculiar feature of ESE systems, apart from complexity, is a memory,
including the long-term one, as well as nonlinear and unstable nature of
interactions and components, which makes it harder to predict their future
behavior.

Unfortunately, mathematical models based on differential equations have no
memory (there is no aftereffect), while for models with memory, where integral
interrelations are used, it is not always possible to take into account
nonlinearity (the integration procedure is linear by definition).

In reality, in the Cauchy problem future systems behavior is defined by its
initial state and doesn’t depend on the way the system reached its current
state. However, it is hardly true that future behavior of a real socio-
economic or socio-ecological system can be predicted by giving an immediate
time “slice” of a variables set that describe its state.

Let us consider possible ways to take into account past events while modeling
ESE systems’ dynamics, which goes beyond the boundaries of classical
differential and integral equations.

Functional differential lagging equation can serve as a simple example of the
dynamic model with memory, where present time is defined by the state variable
$x(t)$ and depends on the past state $x(t-\tau)$ with constant time lag
$\tau=const$:

$x(t)=f\left(x(t-\tau)\right);t\geq t_{0},$ (1)

where $f(x)$ is the known function, with initial conditions being set for the
half-interval $t_{0}-\tau\leq t<t_{0}$ by the function $\phi(t)$:

$x(t)=\phi(t);t_{0}-\tau\leq t<t_{0}.$ (2)

Given the 2 equation 1 has the only solution, defined by recurrent ratios:

$x(t)=\begin{cases}f\left(\phi(t-\tau)\right);\text{ if }t_{0}\leq
t<t+\tau;\\\ f\left(f\left(\phi(t-\tau)\right)\right);\text{ if
}t_{0}+\tau\leq t<t+2\tau;\\\
f\left(\left(f\left(\phi(t-\tau)\right)\right)\right);\text{ if
}t_{0}+2\tau\leq t<t+3\tau;\\\ \cdots\end{cases}$ (3)

Using Dirac delta function, as defined by ratios:

$\delta(t)=0,\text{ if }x\neq
0;\int\limits_{-\infty}^{+\infty}{\delta(t)dt}=1,$ (4)

we can formally rewrite equation 1 in the integral form:

$x(t)=\int\limits_{-\infty}^{t}{dt_{1}f\left(x(t_{1})\right)H(t_{1},t)};H(t_{1},t)\equiv\delta\left(t_{1}-(t-\tau)\right);t\geq
t_{0}.$ (5)

Delta function is not a function in the conventional interpretation and is
related to the class of generalized functions that were mathematically
described only in the middle of the last century r014SobloevFunctProstr
(physics started using this function much earlier). Its classical form is
considered to be a limit of the “peak” sequence, with its centre set in the
point of origin. The afore-mentioned “peaks” indefinitely converge widthway,
indefinitely increase throughout the height and have a unit area.

An approximate classic integral analogue of the equation 5 can be derived by
substituting $\delta(t)$ with an ordinary function - some specific narrow
enough “peak” of a unit area, a certain finite width $\backsim\Delta t$ as
well as a finite height . The derivative of the Fermi function is one of the
possible examples:

$\Phi(t)=\frac{1}{1+exp\left(\frac{-t}{\theta}\right)};\delta(t)\approx\frac{1}{\theta\left(2+exp\left(\frac{-t}{\theta}\right)+exp\left(\frac{t}{\theta}\right)\right)}.$
(6)

If the system’s state in the moment $t$, $x(t)$, is defined not by one, as in
1, but $k$ ($k=2,3,4,\cdots$) of her past states
$x(t-\tau_{1}),x(t-\tau_{2}),\cdots x(t-\tau_{k})$ in the following moments of
time $(t-\tau_{1}),(t-\tau_{2}),\cdots,(t-\tau_{k})$ respectively
$(\tau_{1}=const,\tau_{2}=const,\cdots,\tau_{k}=const,\tau_{1}>\tau_{2}>\cdots>\tau_{k}>0)$,
then instead of (1), (2), (5) we get:

$x(t)=f\left(x(t-\tau_{1});x(t-\tau_{2});\cdots;x(t-\tau_{k})\right);t\geq
t_{0};$ (7) $x(t)=\phi(t);t_{0}-\tau_{1}\leq t<t_{0};$ (8)
$\begin{array}[]{c}x(t)=\int\limits_{-\infty}^{t}{dt_{1}}\int\limits_{-\infty}^{t}{dt_{2}}\cdots\int\limits_{-\infty}^{t}{dt_{k}f\left(x(t_{1}),x(t_{2}),\cdots,x(t_{k})\right)}.\\\
\delta\left((t_{1}-(t-\tau_{1})\right)\delta\left((t_{2}-(t-\tau_{2})\right)\cdots\delta\left((t_{k}-(t-\tau_{k})\right);t\geq
t_{0}.\end{array}$ (9)

Therefore if the system’s state in the time moment $t$ depends on the infinite
sequence of its past states, the integral analogue of the functional
differential lagging equation will, generally speaking, contain an integral of
the infinite multiplicity. At the same time the infinite amount of past states
can relate to both finite $(t-\tau_{1};t)$ (short-term memory) and infinite
$(-\infty;t)$ (long-term memory) time span.

Pay attention that the classic integral lagging equation is one of the
Volterra type r015HeliFunctDiffury :

$x(t)=\int\limits_{-\infty}^{t}{F\left(x(\tilde{t});t;\tilde{t}\right)d\tilde{t}},$
(10)

where $F\left(x(\tilde{t});t;\tilde{t}\right)$ \- is an arbitrary (generally
nonlinear) function of variables $x(\tilde{t});t;\tilde{t}$, which allows to
take into account system’s memory of its past states only in the additive
approximation, which becomes evident, if the right section 10 is rewritten in
the following way:

$\begin{array}[]{c}\int\limits_{-\infty}^{t}{F\left(x(\tilde{t});t;\tilde{t}\right)d\tilde{t}}\equiv\int\limits_{t_{1}}^{t}{F\left(x(\tilde{t});t;\tilde{t}\right)d\tilde{t}}+\int\limits_{t_{2}}^{t_{1}}{F\left(x(\tilde{t});t;\tilde{t}\right)d\tilde{t}}+\cdots=\\\
F\left(x(\tilde{t_{1}});t;\tilde{t_{1}}\right)\cdot(t-t_{1})+F\left(x(\tilde{t_{2}});t;\tilde{t_{2}}\right)\cdot(t_{1}-t_{2})+\cdots;\\\
t>t_{1}>t_{2}>\cdots;\tilde{t_{1}}\in\left[t_{1},t\right];\tilde{t_{2}}\in\left[t_{2},t_{1}\right];\cdots\end{array}$
(11)

In connection with it note that the equation 9 in case of an additive
dependency of contemporaneity on the past, i.e. in case:

$f\left(x(t-\tau_{1});x(t-\tau_{2});\cdots\right)\equiv
f_{1}\left(x(t-\tau_{1})\right);f_{2}\left(x(t-\tau_{2}))\right);\cdots$ (12)

becomes a particular case of the equation 10 with the following integrand:

$F\left(x(\tilde{t});t;\tilde{t}\right)\equiv
f_{1}\left(x(\tilde{t})\right)\delta\left(\tilde{t}-(t-\tau_{1})\right)+f_{2}\left(x(\tilde{t})\right)\delta\left(\tilde{t}-(t-\tau_{2})\right)+\cdots$
(13)

Meaningful analysis of nonlinear models dynamics with memory, in which the
future is defined by the infinite amount of states in the past is generally
possible only in case of a discrete representation. The results of such
analysis will be approximated, i.e. will contain uncertainty, which has to be
considered endogenous, i.e. internal, and peculiar to this very system.

With a certain level of time sampling, models with memory both 7 and 10
becomes:

$x(n+1)=f\left(x(n);x(n-1);x(n-2)\cdots\right).$ (14)

To take into account and quantify the uncertainties, observed in ESE as well
as other complex systems probability models are normally used. However their
application is based on doubtful hypotheses, while the statistical
interpretation of the results is not always informative enough and results
might not correspond with the real process occurring within the system. In
particular, the well-known problem of $1/f$–noise (look for example
r016BukingemShumy ), closely connected to the presence of long-term memory in
complex systems, implies the absence of the mean temporary value (as a limit
of a certain time span converging to infinity, which serves as the basis for
averaging) for any process occurring in such kind of system. Therefore such
processes can’t have a rigorous statistical substantiation.

##  5 Modern concepts in ESE systems modeling

New approaches to modeling and prediction of complex nonlinear systems
dynamics with memory are based on the use of determined chaos and neural
networks technologies (cf. e.g. r011EzhovShumsky98 ; r017LorenzNonlinear89 ;
r018PetersChaosOrder ). Both investigation and realization of such techniques
has become possible only with the appearance of quick-operating computers. Use
of the recurrent computational process has become the general feature for all
these technologies:

$x_{n+1}=f_{n}\left(f_{n-1}\left(\cdots\left(f_{1}(x_{1})\cdots\right)\right)\right),n=1,2,\cdots,$
(15)

where $f_{i}(x_{i})$ is a certain nonlinear mapping of a multi-dimensional
vector $x_{i}$, $i$ \- discrete, real or fictitious, time. Identification of
the model 15 is reduced to the determination of functions $f_{i}(x_{i})$,
while the differences between the models of determined chaos and neural
networks are connected with the function type and methods of its definition
(neural network models normally use rather narrow class of $f_{i}(x_{i})$
mappings r012Zaencev ). Generally speaking, stability or convergence of the
process 15 is not required, whereas a single-step set of vector $x_{i}$
components as well as their time dynamics can be of great interest.

For the particular case of the model 15, introducing corresponding lagged
variables, a model 14 can be transformed.

Both determined stable processes, described by integro-differential equations,
and random processes, which also include complex Markov chains (CMC), can be
formally considered as separate extreme cases of determined chaos models
realization 15. Given the sampling scale, which tends to zero, if such a
tendency makes sense and corresponding limits exist, we derive classical
differential and integral problem statement. Finite $\Delta t$ allows to get
models with discrete time, which in the general case in the corresponding
phase space (which also includes lagged variables) can produce both measurable
sets (discrete or continuous) that allow probabilistic interpretation and
those of the special structure – fractals r019FederFractals , that can’t be
always interpreted in that way.

Various digital generators of so-called random sequences used in imitational
modeling can be an example of determined chaos models that allow probabilistic
interpretation.

Let us note that in reality there are no accurate procedures that would give
an opportunity to distinguish a “real” random sequence from the pseudorandom
one.

##  6 Markov chains prediction technology

Suppose there is a sequence of a certain system discrete states. From this
sequence we can determine transitions probabilities between the two states.
Simple Markov chain is a random process, in which the next state probability
depends solely on the previous state and is independent from the rest of them.
Complex Markov chain, unlike the simple one, stands for the random process, in
which the next state probability depends not only on the current, but also on
the sequence of several previous states (history). The amount of states in
history is the order of the Markov chain.

Theory of simple Markov chains is widely presented in literature, for example
r020TihonovMironovMarkovProcesses . As for the high order Markov chains,
modern literature r021KornVMSpravochnik73  can offer us a mere definition.
Developing complex or high order Markov chain’s properties is not widely
presented in modern scientific publications. It’s necessary to mention the
papers RafteryHighOrderMarkovChain ; RafteryTavare94  where properties of
complex Markov Chains are developed, but no prediction algorithm is proposed
there. The development of prediction method, based on complex Markov chains,
is proposed in this paper.

Markov chain of the higher order can be brought to a simple Markov chain by
introducing the notion of a “generalized state” and including a series of
consequent system’s states into it. In this case, tools of simple Markov
chains can be applied to the complex ones.

Investigated dynamic series is a result of a certain process. It is assumed
that this process is determined, which implies the existence of a causal
dependence of further states on history. It is impossible to fix and analyze
the infinite history, which puts obstacles in the way of an accurate detection
of this influence and making precise predictions.

The problem consists in the maximal use of information, which is contained in
the known segment of the time series, and subsequent modeling of the most
probable future dynamics scenario.

The observed process is described as a time series of prices $p_{t}$ with the
given sampling time span $\Delta t$

$p_{ti}=p(t_{0}+i\Delta t).$ (16)

Discrete presentation of the time series is in fact a way of existence of this
very system. New prices are formed on the basis of contracts or deals, made on
the market in certain discrete moments of time, while the price time series is
a series of the averaged price levels during the chosen time intervals. While
making a decision each trader, who is an active part of the pricing system,
works solely with discrete series of the chosen time interval (e.g. minute,
5-minute, hourly, daily etc.). For $\Delta t\to 0$ the accuracy of data
presentations reaches a certain limit, since for relatively small $\Delta t$
the price leaps in the moment of deal, while staying unchanged and equal to
the last deal during the time between the two deals. Hence, the discreteness
of time series has to be understood not only as a limited presentation of
activity of the complex financial system, but also as one of the principles of
its operation r007KurbanovSaptsin07 ; r008SaptsinMarkovPaper09_ ;
r022SaptsinSoloviev_ ; r022SapSolArxiv_ ; r023SolDerbMonogr2010_ .

The time series of initial conditions has to be turned into a sequence of
discrete states. Let us denote the amount of chosen states as $s$, each of
them being connected to the change in the quantity of the initial signal
(returns). For example, consider the classification with two states, first of
which corresponds to positive returns as the price increases, while the second
one – to negative as it descends. Generally all possible increments of the
initial time series are divided into $s$ groups. Ways of division will be
discussed further.

Next we develop predictions for the time series of sampled states. For the
given order of the Markov chain and the last generalized state the most
probable state is chosen to be the next one. In case if ambiguity occurs while
the state of maximum probability is being evaluated, an algorithm is used that
allows reducing the amount of possible prediction scenarios. Therefore we get
the series of predicted states that can be turned into a sampled sequence of
prognostic values.

Evaluation of increments, prediction, and subsequent restoration are conducted
for the given hierarchy of time increments $t$. To use the given information
as effectively as possible, the prediction is conducted for time increments
$t=1,2,4,8,...$, or a more complex hierarchy of increments and subsequent
“splicing” of the results derived from different prediction samplings.

The procedure of prediction and splicing is iterative and conducted starting
from smaller increments, adding a prediction with the bigger time increment on
every step.

As the sampling time step $t$ increases, the statistics for the investigation
of Markov chains decreases, whereas the biggest sampling step, which takes
part in the prognostication, limits itself. To supplement the prediction with
the low-frequency component the approximation of zero order is being used in
the form of a linear trend or a combination of a linear trend and harmonic
oscillations r024SapChFourierKharkov ; r025ChabS10_ .

##  7 Prediction construction algorithm

Let us consider the consequence of operations, required for the prognostic
time series construction. To do this we need to set the following parameters:

1) The type of time increments hierarchy (simple – powers of two, complex –
product of powers of the first simple numbers).

2) Values of $s$ – the amount of states and $r$ – the order of the Markov
chain. These parameters can be individual for every sampling level; finding of
optimal parameters is done experimentally.

3) Threshold values $\delta$, and minimal number of transitions $N_{min}$.

Prediction construction algorithm includes the following steps:

1) Generating hierarchy of time increments - $t$ sequence. The maximal of them
has to correspond to the length of a prognostic interval $N_{max}$.

2) For every time increment $\Delta t$, as the increments increase, a
prediction of states and restoration of the time series along the prognostic
states is conducted. Current stage includes following actions:

2.1. Evaluating increments (returns) of the series with $\Delta t$ sampling.

2.2. Transforming the time series of increments into the series of state
numbers ($1..s$).

2.3. Calculating transition probabilities for generalized states.

2.4. Constructing the series of prognostic states using the procedure of
defining the most probable next state.

2.5. Restoring the value series from the state series with $\Delta t$
sampling.

2.6. Splicing the prediction of $\Delta t$ sampling with the time series
derived from splicing of the previous layers (with the lesser step $\Delta
t$). In case if the current time series is the first one, the unchanged time
series will come as a result of splicing.

3) To splice the last spliced time series with the continuation of the linear
trend, created along all previously known points.

The time series, spliced with the linear trend, is the result of prediction.
Let us consider the stages of the given algorithm in detail.

##  8 States in complex Markov chains and approaches for defining them

In everything that concerns current technology, states are connected to the
measuring of a prognostic value. There is a number of ways to classify returns
in states, from which the following are suggested. One of them is the
classification based on the homogeneity principle as concerning the amount of
representatives in classes; based on the homogeneity principle of deviation,
as well their combinations for different deviation modules.

Increment or returns of the time series serves as the basis for states
classification r025ChabS10_ ; r026solovievmatheconomics_ . Absolute $r_{a}$
and relative $r_{t}$ increments of the time series are considered:

$r_{a}=p_{t}-p_{t-\Delta t},$ (17) $r_{t}=\frac{p_{t}-p_{t-\Delta t}}{p_{t}},$
(18)

where $p_{t}$ – is the input time series of price dynamics, $\Delta t$ –
sampling interval, which is chosen for subsequent analysis. It is known that
mathematical expectation of the returns time series equals zero, whereas
variation comes as the measure of time series volatility. Based on returns
values $r_{t}$ classification and transformation of values to the time series
of discrete states are conducted. One of the classification principles is
homogeneity according to the amount of class representatives. This
classification divides the set of all increments into $s$ groups equal in
number. Calculated with the given sampling, time series increments are then
systematized in growth and divided into equal parts. Thus we define limit
values $\\{r_{lim,i}\\}$, which are used afterwards during transformation of
the returns into class numbers. Large number of identical states can cause
certain problems, such as identical bounds of several neighbouring states. It
creates a number of states with no representatives, which makes correction of
the division a necessary action. In that way we will reach the largest
possible homogeneity in state division. Classification is conducted along the
following algorithm r027SolSapChDrezden09 ; r028ChabRiga10 :

$s_{t}=(i\mid r_{lim,i-1}>r_{t}>r_{lim,i})$ (19)

where $s_{t}$ is the number of state, which corresponds to the moment of time
$t$, for which the returns level was computed $r_{t}$; i is the number of
state $[1\dots s]$, which is characterized by the interval
$[r_{lim,i-1},r_{lim,i}]$ corresponding to the calculated returns level
$r_{t}$.

Apart from the returns interval, given by the aforementioned values
$[r_{lim,i-1},r_{lim,i}]$, a mean returns value is chosen for every state
$r_{avg,i}$, which will be used in time series values transformation according
to the prognostic discrete states.

Another way of dividing the time series into states implies dividing the
interval of returns values into equal parts, from minimal to maximal
deviation. In this case homogeneity according to the amount of representatives
in states does not occur. In fact this method differs from the previous one in
terms of defining limit values $\\{r_{lim,i}\\}$. Possible combined ways of
division, in case of which the limit value, dependent on standard deviation,
is used instead of maximal and minimal value, and division is conducted
homogenously according to the deviation.

Since the real causal dependence is unknown during the process, to find it
adequate state classification, which would allow to reveal vital dependencies
of the time series, is required. We suggest a couple of ways to divide the
time series into states, which in the first place allow to preserve adequate
transition probabilities between states, as well as prevent averaged
deviations inside the states from affecting the accuracy of the derived
prediction.

To check the efficiency of division we conduct the sampling procedure and
classify the increments according to each hierarchy. Having completed that, we
restore the time series using known states for each hierarchy and finish the
splicing procedure. Since the state series correspond to the initial time
series, we get the curve, with deviation, caused exclusively by the state
averaging mistake (quantum mistake). Thus, having set a certain value of state
numbers $s$ and carried out sampling, restoring, and sampling procedures
(excluding prediction), we get absolute sampling (quantum) mistake.

Increasing the number of states, we improve the accuracy of restoration,
however one should remember, that the choice of the quantum levels is limited
by the fact that the transition probabilities definition with sufficient
accuracy is required, which is confirmed by artificial test time series
prediction experiments.

##  9 Step-by-step prediction procedure. Defining the most probable state on
the next step, prediction scenarios

Predicting procedure uses the most probable state as the next one under
current circumstances. Probability matrix of state transitions is used for the
afore-mentioned purpose. In this case, you have to take into account that
probabilities are calculated with a certain mistake. We cannot precisely
compute the probabilities, since it is impossible to derive an infinite time
series, and only a part of the time series is known – the known part serves as
the basis for probabilities. The second important aspect implies the case of
several states with maximal probability.

To prevent the omission of the states, for which the probabilities are
computed with a mistake, one should add a state with maximal probability to
the states, which are located in the distance of $\delta$ from the maximal
one. The value of parameter $\delta$ depends on the probability evaluation
mistake and requires experimental refinement.

If $\delta>0$, the number of states with maximal probability increases in
comparison to the value $\delta=0$. Let us call a couple of neighbouring
states with maximal probability a cluster. Cluster states with average
deviation values are supposed to have the largest probability. To predict the
dynamics, let us confine ourselves to one or two most probable states. To
define them a following algorithm is suggested:

1) If levels (discretized increments) create several clusters (cluster is a
group of several neighbouring levels – cluster elements, minimal cluster is a
single isolated level) with maximal probability, we choose the largest
cluster.

2) If the number of cluster elements is odd, as $k_{max}$ we choose a central
element.

3) If the number of cluster elements is even, we consider two central cluster
elements and choose as $k_{max}$ the one, which is closer to the centre of
distribution.

4) If two central cluster elements are equidistant to the centre of
distribution, we consider both cases as possible variants of $k_{max}$ values
(bifurcation point).

5) If there are several clusters of maximal size, we consider them as new
elements, which can also form clusters that will undergo the same steps 1)-4).

This principle is based on the following ideas:

1) If there are two neighbouring states of maximal probability, it is better
to take the one, which is closer to the centre of distribution, in order to
minimize the risk of occurrence of false linear trends in the prediction.

2) If levels of maximal probability are not the neighbouring ones, at least
two variants have to be considered, as it can be connected to the bifurcations
that should not be omitted.

3) If the prediction is carried out according to 1) (on all stages of the
hierarchy), we receive a certain approximation of the lower limit of the
prediction, whereas in case of 2) – we get an approximation of the upper
limit.

Hence this algorithm can adequately restore the case of possible bimodal
probability distribution, it is proposed to consider 2 prediction scenarios.

In case of the complex Markov chains, probability of the next state depends
not only on the previous state, but also on the sequence of $r$ states, which
have occurred before given. In this case, it is necessary to calculate
transition probabilities from the sequence of $r$ states into the $r+1$ state.
Formally, these probabilities can be written into the rectangular table of
$(r^{s},s)$ size.

Having generalized the notion of “present state” and included a sequence of
$r$ preceding states into it, we can reduce Markov chains of $r$ order to the
chain of the first order. Thus transition probabilities can be written into
rectangular matrices of $(r^{s},r^{s})$, that come as transition probability
matrices for generalized states.

The process of prediction implies the following: the last state is chosen (in
case of Markov chains of an order $r>1$ a sequence of $r$ latest states is
taken). The probability of transition from current state to all possible
states is defined. From all possible states a state with maximal probability
is chosen. It is possible that several states with maximal probability occur,
which can be explained by the bimodal probability distribution. The process of
decision-making in this case is described later.

The chosen most probable state is taken as the next prognostic state and the
procedure is repeated for the next (last added) state. Thus we receive a time
series of prognostic states for the given sampling time $\Delta t$.

Further according to the received state sequence and known initial value the
time series is being restored for the given time sampling $\Delta t$. In this
case every state implies $\Delta t$ points of the time series. On the stage of
state classification every state was connected to the average increment
$r_{avg,i}$, which is added to the value of the last point in the time series,
and the next discrete point is computed. Intermediate points are filled as
linear interpolation of two known neighbouring points. Algorithm of $y_{t}$
time series values restoration according to the initial price $p_{t}$ and a
series of average increments $r_{avg,ik}$, corresponding to the prognostic
states $s_{k}$, can be given by a sequence of calculations:

$\begin{array}[]{c}y_{t}=p_{t},\\\ y_{t+1}=y_{t}+r_{avg,i1}/\Delta
t=p_{t}+r_{avg,i1}/\Delta t,\\\ y_{t+2}=y_{t+1}+r_{avg,i1}/\Delta
t=p_{t}+2r_{avg,i1}/\Delta t,\\\ \dots\\\ y_{t+\Delta t-1}=y_{t+\Delta
t-2}+r_{avg,i1}/\Delta t=p_{t}+(\Delta t-1)r_{avg,i}/t,\\\ y_{t+\Delta
t}=y_{t+\Delta t-1}+r_{avg,i1}/\Delta t=p_{t}+\Delta tr_{avg,i}/\Delta
t=p_{t}+r_{avg,i1},\\\ y_{t+\Delta t+1}=y_{t+\Delta t}+r_{avg,i2}/\Delta
t=p_{t}+r_{avg,i1}+r_{avg,i2}/\Delta t,\\\ \dots\\\ y_{t+n\Delta
t-1}=y_{t+n\Delta t-2}+r_{avg,in}/\Delta
t=p_{t}+\sum_{k=1}^{n-1}r_{avg,ik}+\frac{\left(\Delta t-1\right)}{\Delta
t}r_{avg,ik},\\\ y_{t+n\Delta t}=y_{t+n\Delta t-1}+r_{avg,in}/\Delta
t=p_{t}+\sum_{k=1}^{n}r_{avg,ik}.\\\ \end{array}$ (20)

##  10 Time increments hierarchy and splicing procedure

Time series increments will be computed with different steps. For example,
analogous to the discrete Fourier transform, time increments are equal the
powers of 2 are considered. First, we calculate increments as a remainder of
two nearest neighbouring time series values, then next nearest values are
considered with the step of 2, 4, 8, 16 etc. Let us mark this difference in
time as $\Delta t$.

For every $\Delta t$ we conduct an increment time series transformation
leading to a time series of states. Further we predict the future sequence of
states and restore the time series with the given sampling rate according to
the prognostic series of states.

Time series, received as a result of restoring for different $\Delta t$,
undergo the splicing procedure, which gives out an actual prognostic time
series.

Thus an increment hierarchy is chosen, where each one is responsible for its
own sampling rate, which serves as a basis for predicting, restoring and
splicing.

The splicing process implies the following. The procedure is iterative. With
every next (along with the increasing step) sampling time the series corrects
itself, driving the prediction, formed under lower $\Delta t$, to its actual
point. Transformations that are conducted during splicing can be written down
in the form of the following calculations.

Suppose the splicing procedure has been finished for all time increments
$\Delta t<Deltat_{i}$, the prediction has been done under the $\Delta t_{i}$
sampling according to formulae 20, and as a result a time series $y_{i}$ has
been derived. Let us consider the iterative splicing procedure of the received
series $y_{i}$ with the series, acquired during all preceding splicing
procedures $g_{i}$.

Since the series $y_{i}$ contains system points only in moments aliquot to
$\Delta t_{i}$, and other points of the series are interpolated, the process
of splicing implies the substitution of these interpolated points with the
values of system points from previous $\Delta t<\Delta t_{i}$, which are
contained in the series of results of previous splicing procedures $g_{i}$.
Splicing algorithm can be written in the sequence of computations:

$\begin{array}[]{c}z_{t}=g_{t}=p_{t},\\\ z_{t+1}=g_{t+1}+\left(y_{t+\Delta
t_{i}}-g_{t+\Delta ti}\right)/\Delta t_{i},\\\
z_{t+2}=g_{t+2}+2\left(y_{t+\Delta t_{i}}-g_{t+\Delta ti}\right)/\Delta
t_{i},\\\ \dots\\\ z_{t+\Delta t-1}=g_{t+\Delta t-1}+\left(\Delta
t_{i}-1\right)\left(y_{t+\Delta t_{i}}-g_{t+\Delta ti}\right)/\Delta t_{i},\\\
z_{t+\Delta t}=g_{t+\Delta t}+\left(\Delta t_{i}\right)\left(y_{t+\Delta
t_{i}}-g_{t+\Delta ti}\right)/\Delta t_{i}=y_{t+\Delta_{ti}},\\\ z_{t+\Delta
t+1}=g_{t+\Delta t+1}+\left(\left(y_{t+2\Delta ti}-g_{t+2\Delta
ti}\right)-\left(y_{t+\Delta ti}-g_{t+\Delta ti}\right)\right)/\Delta
t_{i},\\\ z_{t+\Delta t+2}=g_{t+\Delta t+2}+2\left(\left(y_{t+2\Delta
ti}-g_{t+2\Delta ti}\right)-\left(y_{t+\Delta ti}-g_{t+\Delta
ti}\right)\right)/\Delta t_{i},\\\ \dots\\\ z_{t+n\Delta t-1}=g_{t+n\Delta
t-1}+\frac{\left(\Delta t-1\right)}{\Delta t}\left(\left(y_{t-n\Delta
t}-g_{t-n\Delta t}\right)-\left(y_{t-(n-1)\Delta t}-g_{t-(n-1)\Delta
t}\right)\right),\\\ z_{t+n\Delta t}=g_{t+n\Delta t}+\left(\left(y_{t-n\Delta
t}-g_{t+n\Delta t}\right)-\left(y_{t+(n-1)\Delta t}-g_{t+(n-1)\Delta
t}\right)\right)=\\\ =g_{t+(n-1)\Delta t}-y_{t+(n-1)\Delta t}-y_{t-n\Delta
t}.\end{array}$ (21)

##  11 Results of stock indices prediction

In this section we offer the results of stock indices prediction. The stock’s
indices databases are available from Finance_yahoo . Point 2000 indicates the
starting moment of the prognosis: March 24, 2011. The green line on the next
figures indicates real indice’s or price’s values. Our software for time
series forecasting by the proposed methods is available from our website:
http://kafek.at.ua/MarkovChains1_2_20100505.rar.

Figure 1: Stock indices prediction. a) Dow Jones Industrial Average - DJI
(USA). b) FTSE 100 (Great Britain)

Figure 2: Financial companie’s share prices forecasting. a) Morgan Stanley
(USA). b) BNP Paribas (France)

Prediction time series with different input learning set’s length are shown at
the fig.3 and 4. Prediction series for DJI at the figure 3 are more
correlated, than FTSE index at the figure 4. At the subplot b) of the above
mentioned plots the mean value and standard deviations of the prediction’s
series are presented. The time of prediction series beginning on the next
figures is the point 1000 and correspond to October 14, 2011.

Figure 3: Dow Jones Industrial Average - DJI (USA). a) Prediction series,
calculated with different learning set’s length. b) Mean value and standard
deviation for prediction series.

Figure 4: FTSE 100 index prediction. a) Prediction series, calculated with
different learning set’s length. b) Mean value and standard deviation for
prediction series.

The normalization procedure is proposed in order to compare indices and it’s
prediction series with different absolute values. The normalized values
calculated with the following formula:

$y_{n}(t)=\frac{y(t)-min\left(y(t)\right)}{max\left(y(t)\right)-min\left(y(t)\right)}.$
(22)

Normalized prediction time series are shown at the fig.5 (America), fig.7
(Europe, developed countries), fig.7 (Europe, PIIGS), fig.9 (Asian markets).
All the figures contain mean time series, which are weighted average of
countrie’s stock indices predictions, weihted with GDP values EconomyWatch_com
for the corresponding countries.

Figure 5:  Normalized mean values for the prediction series of America’s stock
indices. Brazil (BVSP), Mexico (MXX), Canada (GSPTSE), Argentina (MERV), USA
(S&P 500)

Figure 6:  Normalized mean values for the prediction series of European stock
indices. Developed countries: FTSE (Great Britain), DAX (Germany) FCHI
(France), Netherlands (AEX)

Figure 7:  Normalized mean values for the prediction series of European stock
indices. Portugal (PSI20), Italy (FTSEMIB), Ireland (ISEQ), Greece (GD) and
Spain (IBEX).

Figure 8:  Normalized mean values for the prediction series of Asian stock
indices. China (SSEC, HSI), Korea (KS11), Japan (NIKKEI), India (BSESN), New
Zealand (NZ50).

Figure 9:  Mean values of normalized World’s powerful economies indices
prediction series.

##  12 Conclusions and further work

Current paper suggests an algorithm of time series prediction based on complex
Markov chains. Hierarchy of time increments principle allows to use the
information, which is contained in the time series during the prognosis
construction, to its fullest. Experimental work on stock market indices time
series prediction shows the efficiency of the algorithm and confirms the
relevance of further research of the offered method.

## References

  * [1] A. A. Samarskii and A. P. Mikhailov.  Mathematical Modelling: Ideas. Methods. Examples.  Fizmatlit, Moscow, 2001.
  * [2] V. V. Bohoboyaschyy, K. R. Kurbanov, P.B. Paly, and V.M. Shmandiy.  Principles of Forecasting in Ecology: Textbook (in Ukrainian).  Center navchalnoji literaturi, Kyiv, 2004.
  * [3] O. G. Ivakhnenko.  Grouping method of data handling - the concurrent of stochastic approximation methods (in ukrainian)  Automatika, 3(3):58–72, 1968.
  * [4] J. V. Andersen, S. Gluzman, and D. Sornette.  Fundamental framework for “technical analysis” of market prices.  European Physical Journal B, 14:579–601, March 2000.
  * [5] D. Horvath and R. Pincak.  From the currency rate quotations onto strings and brane world scenarios. arXiv:1104.4716 , April 2011.
  * [6] M. Repasan and R. Pincak.  The string prediction models as application to financial forex market. arXiv:1109.0435 , September 2011.
  * [7] P. V. Elyutin and V. D. Krivchenkov.  Quantum Mechanics with tasks.  Nauka, Moscow, 1976.
  * [8] L.D. Landau and E.M. Lifshitz.  Quantum mechanics: non-relativistic theory.  Teoreticheska fizika (Izd. 3-e) (Landau, L. D, 1908-1968). Butterworth-Heinemann, 1977.
  * [9] J. Sapir.  K Ekonomitcheskoj teorii neodnorodnyh sistem - opyt issledovanija decentralizovannoj ekonomiki.  GU VShE, Moskov, 2001.
  * [10] L. von Bertalanffy.  General system theory—a critical review.  “General Systems”, VII:1–20, 1962.
  * [11] K. R. Kurbanov and V. M. Saptsin.  Markov chains as technology for social, economic and ecological processes forecasting.  In “Problemy ta perspectivy rozvitku reg onalno rinkovo ekonom ki” Conference proceedings (in Russian), pages 10 – 14, Kremenchuk, May, 11-13 2007.
  * [12] V. M. Saptsin.  Experience of using genetically complex Markov chains for the neural network technology forecasting.  Visnyk Krivorizkogo ekonomichnogo institutu KNEU, 2 (18):56 – 66, 2009.
  * [13] Y. P. Lukashin.  Adaptive Methods of Time Series Forecasting: Textbook.  Finance and Statistics, Moscow, 2003.
  * [14] Y. P. Zaichenko.  Fuzzy models and techniques in intelligent systems. Monograph(in Russian).  Slovo, Kiev, 2008.
  * [15] A. A. Ezhov and S. A. Shumsky.  Neurocomputing and its application in economics and business (Series “Textbooks” of Economic-Analytical Institute MEPI) / Ed. prof. V. V. Kharitonov.  MEPI, Moscow, 1998.
  * [16] I. V. Zayencev.  Neural networks: basic models.  Textbook for the course “Neural networks” for 5-th grade students (in Russian) Physical Electronics Department, Faculty of Voronezh State University  Voronezh State University, Voronezh.
  * [17] D. M. Chabanenko.  Detection of short- and long-term memory and time series prediction methods of complex Markov chains.  In Visnyk Natsionalnogo tehnichnogo universitetu “Kharkivsky politehnichny institut”. Zbirnik Naukovyh pratz. Tematichny vypusk: Informatika i modelyuvannya (in Ukrainian), number 31, pages 184 – 190. NTU KHPI, Kharkov, 2010. http://www.pim.net.ua/ARCH_F/V_pim_10.pdf
  * [18] S. L. Sobolev.  Selected topics from the theory of functional spaces and generalized functions.  Nauka, Moscow, 1989.
  * [19] J. Hale.  Theory of Functional Differential Equations.  Springer-Verlag, New York, Heidelberg, Berlin, 1977.
  * [20] M.J. Buckingham.  Noise in electronic devices and systems.  Ellis Horwood series in electrical and electronic engineering. E. Horwood, 1983.
  * [21] Hanz-Valter Lorenz.  Nonlinear Dynamical Economics and Chaotic Motion.  Springer-Verlag, 1989.
  * [22] E.E. Peters.  Chaos and order in the capital markets: a new view of cycles, prices, and market volatility.  Wiley finance editions. Wiley, 1996.
  * [23] J. Feder.  Fractals.  Physics of solids and liquids. Plenum Press, 1988.
  * [24] V. I. Tikhonov and V. A. Mironov.  Markov Processes.  Soviet Radio, Moscow, 1977.
  * [25] G.A. Korn and T.M. Korn.  Mathematical handbook for scientists and engineers: definitions, theorems, and formulas for reference and review.  Dover books on mathematics. Dover Publications, 2000.
  * [26] Andrian E. Raftery.  A model for high-order markov chains.  Journal of the Royal Statistical Society., 1985.
  * [27] Adrian Raftery and Simon Tavare.  Estimation and modelling repeated patterns in high order markov chains with the mixture transition distribution model.  Appl. Statist., 43(1):179–199, 1994.
  * [28] V. M. Saptsin and V. N. Soloviev.  Relativistic Quantum Econophysics. New paradigms of Complex systems modeling: Monograph. http://kafek.at.ua/sol_sap_monogr.rar Brama-Ukraine, Cherkassy, 2009.
  * [29] V. Saptsin and V. Soloviev. Relativistic quantum econophysics - new paradigms in complex systems modelling. arXiv:0907.1142v1 [physics.soc-ph].
  * [30] V. D. Derbentsev, A. A. Serdyuk, V. N. Soloviev, and O. D. Sharapov.  Synergetical and econophysical methods for the modeling of dynamic and structural characteristics of economic systems. Monograph (in Ukrainian).  Brama-Ukraine, Cherkassy, 2010. http://kafek.at.ua/Monogr.pdf
  * [31] V. M. Saptsin and D. N. Chabanenko.  Fourier-based forecasting of low-frequency components of economical dynamic’s time series.  In Problemy ekonomichnoyi kibernetiki: Tezy dopovidey XIV Vseukrayinskoyi Naukovo-praktichnoyi konferentsiyi.(in Ukrainian), page 132, Kharkiv, Oct 8-9, 2009 2009. KhNU imeni VN Karazina.
  * [32] D. M. Chabanenko.  Discrete Fourier-based forecasting of time series.  Sistemni tehnologii. Regionalny mizhvuzivsky zbirnik naukovyh pratz (in Ukrainian), 1(66):114 – 121, 2010. http://www.nbuv.gov.ua/portal/natural/syte/2010_1/15.pdf
  * [33] V. N. Soloviev.  Mathematical economics. A textbook for self-study (in Ukrainian).  CHNU, Cherkassy, 2008.  http://kafek.at.ua/Posibnyk_Soloviev.rar
  * [34] Vladimir Soloviev, Vladimir Saptsin, and Dmitry Chabanenko.  Prediction of financial time series with the technology of high-order markov chains.  In Working Group on Physics of Socio-economic Systems (AGSOE), Drezden, 2009. http://www.dpg-verhandlungen.de/2009/dresden/agsoe.pdf
  * [35] V. Soloviev, V. Saptsin, and D. Chabanenko.  Financial time series prediction with the technology of complex markov chains.  Computer Modelling and New Technologies, 14(3):63–67, 2010. http://www.tsi.lv/RSR/vol14_3/14_3-7.pdf
  * [36] Yahoo! finance database. http://finance.yahoo.com
  * [37] Economic indicators : Gdp per capita (current prices, national currency) 2010. http://www.economywatch.com/economic-statistics/economic-indicators/GDP_Per_Capita_Current_Prices_National_Currency/2010/


# Band gap tunning in BN-doped graphene systems with high mobility

T. P. Kaloni1, R. P. Joshi2, N. P. Adhikari2, and U. Schwingenschlögl2
<EMAIL_ADDRESS>1KAUST, PSE Division, Thuwal 23955-6900,
Kingdom of Saudi Arabia  2Central Department of Physics, Tribhuvan University,
Kirtipur, Kathmandu, Nepal

###### Abstract

Using density functional theory, we present a comparative study of the
electronic properties of BN-doped graphene monolayer, bilayer, trilayer, and
multilayer systems. In addition, we address a superlattice of pristine and BN-
doped graphene. Five doping concentrations between 12.5% and 75% are
considered, for which we obtain band gaps from 0.02 eV to 2.43 eV. We show
that the effective mass varies between 0.007 and 0.209 free electron masses,
resembling a high mobility of the charge carriers.

The outstanding electronic properties of graphene have caused huge interest in
this material, though the gapless features restrict to build graphene-based
electronic devices, which could be alternative for existing silicon-based
technology. Therefore, opening a finite and tunable band gap in graphene is
demanded. Great efforts have been devoted to opening of a tunable gap in
graphene. Various routes have been proposed, for example: superlattices of
graphene with $h$-BN Zhou ; kaloni ; Giovannetti ; Xu , atomic layers of
hybridized $h$-BN and graphene domains ajayan , bilayer graphene park ,
graphene nanowires son , and twisted graphene li . In gated bilayer graphene,
a band gap of 250 meV and in graphene nanowires a band gap of 24 meV has been
opened. The substrate induced band gap is found to be 52 meV. In twisted
graphene a band gap 270 meV have been reported. The many other ways also
proposed to open a band gap in graphene, such as oxidation of graphene Bagri ;
Dai ; Nourbakhsh , hydrogenation Pumera , and inducing of fluorine/boron-oxide
in graphene Withers ; udo2 ; udo3 . Recently, experimentally, a novel and
promising way to open a band gap in graphene is proposed by controlling the
domain size of in-plane heterostructures of graphene and $h$-BN ajayan1

The $h$-BN is wide band gap insulator with the same honeycomb structure as
graphene Watanabe ; Kubota . It exhibits many attractive properties such as
high in-plane mechanical strength, excellent chemical inertness, and high
thermal conductivity Balandin ; Frank ; Song ; Kho ; Chen . It is suitable
material to form the hybrid structure with graphene due to minimal lattice
mismatch of about 1.6% and both having same structural dimensionality. Hence,
it is easy to achieve the hybrid structures of graphene and $h$-BN either in-
plane manner or in stacking scheme as it has been demonstrated experimentally
in Refs. ajayan ; ajayan1 . Experimentally and theoretically the B/N doped
graphene have been studied Sharifi ; Han ; Xue ; sugata  and BN doped graphene
leads to semiconducting Kim ; Leenaerts ; sugata1 . Still mechanism of the
opening of the band gap in BN doped graphene, its electron effective mass, and
mobility are unsettled.

Recently, the BN doped multilayer graphene and BN domain in graphene has been
synthesized Kim ; ajayan1 . Motivated by these experimental reports, we
investigate the effect BN doping in the electronic properties of graphene
monolayer, bilayer, trilayer, multilayer, and superlattice of pristine and BN-
doped graphene, using density functional theory. The opening of the band gap
found to be 0.02 eV to 2.43 eV for 12.5% to 75% BN concentration. The
calculated value of band gap are in good agreement with the available
experimental and theoretical data. We estimate the electron effective mass for
all the systems under study and found to be the value varies form 0.007 me
(which is smaller than that of reported for similar systems) to 0.209 me. We
also estimate the mobility, which is as high as graphene on $h$-BN.

The results presented here are obtained from first-principles density
functional theory in generalized gradient approximation Perdew . The QUANTUM-
ESPRESSO code is used with a plane-wave cutoff energy of 45 Ry paolo . In
addition, Van der Waals interaction has been taken into account via Grimmes
scheme grimme , which is expected to provide correct interlayer spacing and
dispersion. A MonkhorstPack with a 8$\times$8$\times$1 k-mesh is used for all
calculations under study. A 4$\times$4 supercell of graphene is used in our
calculations. For all the structures, the atomic positions are fully relaxed
until the energy convergence of 10-5 Ry and the force convergence of 10-4
Bohr/Ry is achieved.

The BN-doped graphene is classifies as (1) monolayer (see Fig. 1), (2)
bilayer, (3) trilayer, (4) multilayer, and (5) superlattice. In all these
systems, we have used five different BN-doping concentrations such as 12.5%
(2B2N), 25% (4B4N), 43.75% (7B7N), 56.25% (9B9N), and 75% (12B12N). In all of
the mentioned doping concentration the ratio of B to N atom is same ajayan .
Turning to stacked structure: bilayer, trilayer, multilayer (superlattice of
BN-doped graphene), and superlattice of pristine graphene with BN-doped
graphene, there are two stacking schemes, such as AA and AB with different
geometry. All of the possible configurations, N centered orientation in AB
stacking is the most stable Giovannetti . Hence, we only consider N centered
AB stacking structures in our calculations. The distance between two layer in
BN-doped bilayer and trilayer graphene is found to be 3.34 Å which is same as
the separation between two layer in graphite. In all system under study, we
found $a=b=9.94$ Å. We use a vacuum layer of 12 Å in order to avoid the
interlayer interactions due to the periodic boundary conditions. Multilayer
and superlattice of pristine and BN-doped graphene have orientation similar to
that of BN-doped bilayer graphene in which the distance between two layer is
equal to 3.36 Å and the length of the $c$ lattice to be 6.72 Å. The modified
C$-$C bond length for all the systems under study are vary from 1.40 Å to 1.43
Å, C$-$B bond length varies from 1.46 Å to 1.49 Å, C$-$N bond length varies
from 1.36 Å to 1.40 Å, and B$-$N bond length varies from 1.40 Åto 1.45 Å,
while for pristine graphene C$-$C bond length is 1.42 Å and B$-$N bond length
is 1.44 Å for pristine $h$-BN. The bond angle in BN-doped systems varies from
117∘ to 122∘, which is 120∘ for pristine graphene/$h$-BN.

Figure 1: Supercells used for modeling different BN concentrations.

To understand the stability of the structure and to compare the energetics we
calculate the cohesive energy. The cohesive energy is calculated as
$E_{coh}=\frac{E_{cell}-n\cdot E_{C}-m\cdot E_{B}-p\cdot E_{N}}{n+m+p}$, where
Ecoh, cohesive energy, Ecell is the total energy of BN doped 4$\times$4
supercell, EC, EB, and $E_{N}$ are the total energy of isolated C, B, and N
atoms, which is calculated by placing a single atom in a large and n, m, p are
the total number of C, B, and N atoms, respectively. The calculated values of
cohesive energy for all the systems under study are presented in Table I. The
cohesive energy for BN-doped monolayer ranges from 7.76 eV/atom to 8.02
eV/atom for 12.5% to 75%, which value is 7.75 eV/atom to 8.00 eV/atom for both
bilayer and trilayer, and 7.75 eV/atom to 8.12 eV/atom for multilayer and
superlattice of pristine and BN-doped graphene. These value slightly higher
than that of the cohesive energy of 7.41 eV/atom for graphite Greenwood ,
close to superlattices of graphene and $h$-BN kaloni , and also well agree
with the obtained value of the cohesive energy of $h$-BN ciraci . From these
obtained data for cohesive energy, we can conclude that all the system under
study are highly feasible.

The band gap is an intrinsic property of semiconductors, which indeed hugely
determines the transport and the optical properties. It should have a key role
in modern device industry and technology. In this contrast, we focus the
sensitivity of the electronic properties of the BN-doped thicker and thinner
graphene by varying the BN concentration of 12.5% to 75%. The calculated
electronic structures of pristine graphene and BN-doped thinner and thicker
graphene along the $\Gamma$-K-M-$\Gamma$ path is addressed in Fig. 2. The top
row represents a monolayer graphene doped by 12.5% to 75% of BN atoms, second
row represents bilayer, third row represents trilayer, fourth row represents
multilayer, and fifth row represents superlattice of pristine graphene and BN-
doped graphene with the same concentration as monolayer. We obtain the band
gap of 0.42 to 2.43 eV for 12.5% to 75%, respectively, for all the systems
under study, see Table I. The increasing of the band gap with respect to the
BN concentration in graphene is addressed in Fig. 3. The opening of the band
gap is attributed by the broken sublattice symmetry of the pristine graphene
due the presence of the BN pair. Recently, the experimental setup for such
level of doping is established by mass production of graphene samples using
CH4 gas as C source and ammonia borane precursor for BN source, which are
introduced as a vapor in the whole duration of graphene growth Chang . The
concentration of the BN in graphene sample can be varied and controlled by
controlling annealing temperature of the precursor. This process can be
monitored by the X-ray photoelectron spectroscopy, which indeed confirms the
monotonus increment in the BN concentration by the rising of the annealing
temperature of the precursor. The authors of Ref. Chang  measured the band gap
of about 0.6 eV for low concentration of the BN in the graphene sample, which
in fact very close to our calculated band gap for the low concentration of
12.5%. The another reason for the opening of the band gap is the change in the
on-site energy of graphene due to the isoelectronic co-doping of BN pairs Liu1
. The obtained value of the band gap is increasing monotonically with
increasing the BN pairs in the graphene sheet, due to the changing the on-site
energy of the more and more carbon atoms in the system.

The electronic band structure of undoped and the BN-doped bilayer and trilayer
graphene by 12.5% to 75% are shown in second and third rows of the Fig. 2. The
pristine bilayer/trilayer graphene is metallic in nature Zhang1 ; Freitag .
The variation of the band gap follows the similar trend as in the case of
monolayer but with slightly smaller gap, see Table I and Fig. 3. The band gap
is decreasing with increasing the number of layers due the fact that the
inversion symmetry broken for the bilayer/trilayer BN-doped graphene. The band
gap opening in pristine bilayer graphene by the breaking of the inversion
symmetry due to the application of the finite electric field has been
demonstrated experimentally Zhang1 , a band gap of 250 meV has been achieved.
However, in our calculations, by BN-doping in bilayer/trilayer graphene a
variable band gap of 0.04 eV to 2.20 eV is obtained, which can be enhanced and
controlled by applying the electric field as it has been investigated for
pristine bilayer graphene. The system we have studied here have been
synthesized experimentally ajayan ; ajayan1 ; Kim .

Figure 2: Electronic band structure of BN-doped graphene for concentrations
between 0% (red lines) and 75% (blue lines). First row: monolayer, second row:
bilayer, third row: trilayer, fourth row: multilayer, fifth row: superlattice
of pristine and BN-doped layers. Note that the high symmetric points of the
$4\times 4$ supercell are folded back to the same points of the $1\times 1$
cell.

The quantization of the energy states in the potential wells are important to
understand the characteristic features of the superlattices. Which is always
depends on the barrier height, when the barrier hight is small then 2 Å the
electronic energy bands should split into quasi-two-dimansional subbands like
in conventional superlattices, for example InAs/GaSb. Moreover, if the barrier
height is greater than 2 Å, the possibility of the overlapping of the
neighbouring wavefunctions become suppressed. This give rise the the separate
quantum wells of the superlattices Xu . In our case the superlattices having a
barrier height of greater than 2 Å (i.e. 3.36 Å). Our calculated band gap for
superlattice is smaller than that of BN doped monolayer, bilayer, and trilayer
graphene, see the electronic band structure addressed in fourth row of the
Fig. 2 and Table I. Following the same trend as the other cases discussed
above, the band gap increasing with the increasing the concentration of the BN
atoms, which is well agree with the available theoretical prediction Lin  and
experimental observation Kim . Essentially, the electronic band structure of
the superlattice is presented in fifth row of the Fig. 2 and obtained data are
summarized in Table I. In this case, the minimum band gap of 0.08 eV to the
maximum of 0.13 eV is observed. These obtained values are well agree with the
available data for similar systems kaloni ; Giovannetti ; Xu . The
superlattice with all the BN doping concentration would be a promising
candidate to build electronic device in a low voltage limit Zhang1 . All the
systems with low BN doping concentration are good for thermoelectric
applications.

Figure 3: Band gap as a function of doping for various BN-doped graphene
systems.

From tight-binding approximation, the energy spectrum of gapped graphene,
where the charge carriers behaves as massive Dirac particles is expressed as
$E(k)=\pm v\sqrt{k_{x}^{2}+k_{y}^{2}+m^{*}{{}^{2}}v^{2}}$ Casolo , where $\pm$
represents the conduction and valance bands, kx and ky are wave vector and m∗
is electron effective mass and hence band gap can be formulated as
$E_{g}=2m^{*}v^{2}$. Finally, we attempt to calculate the electron effective
mass of the electron defined by $(m^{*})\sim\frac{E_{g}}{2v^{2}}$ apl , where
Eg is the band gap and v is the Fermi velocity, which is found to be
$\sim$0.8$\cdot 10^{6}$ m/s for BN doped systems. It is reported that the
interlayer spacing is the key point to control the band gap in bilayer
graphene but with the increasing the band gap the value of m∗ is hugely
increased, which is not good for device application because it suppress the
carrier mobility. Hence, its very important to keep and control the m∗ as much
as small. Our calculated electron effective mass for all the systems under
study is varies from 0.007 me to 0.209 me, presented in Table I. The
calculated m∗ really very small for low concentration limit of all the BN
doped graphene systems. Interestingly the m∗ is lowest form all the
concentration from 12.5% to 75% for superlattice of pristine and BN-doped
graphene because of smallest band gap opening in this system, and hence would
be great candidate to built a electronic and opto-electronic devices, those
can be operated in a low voltage/energy gap with a highly enhanced mobility
magno .

Quantity | Concentration (%) | Monolayer | Bilayer | Trilayer | Multilayer | Superlattice
---|---|---|---|---|---|---
Ecoh (eV) | 12.5 | 8.02 | 8.00 | 8.00 | 8.12 | 8.18
25 | 7.94 | 7.93 | 7.93 | 7.92 | 8.15
43.75 | 7.85 | 7.84 | 7.84 | 7.83 | 8.11
56.25 | 7.78 | 7.79 | 7.77 | 7.77 | 8.08
75 | 7.76 | 7.75 | 7.75 | 7.75 | 8.07
Eg (eV) | 12.5 | 0.42 | 0.11 | 0.04 | 0.02 | 0.08
25 | 0.85 | 0.50 | 0.40 | 0.49 | 0.12
43.75 | 1.46 | 1.30 | 1.15 | 0.96 | 0.20
56.25 | 1.76 | 1.52 | 1.40 | 1.22 | 0.14
75 | 2.43 | 2.20 | 2.08 | 1.90 | 0.13
m∗ (me) | 12.5 | 0.040 | 0.009 | 0.003 | 0.02 | 0.007
25 | 0.073 | 0.043 | 0.034 | 0.042 | 0.010
43.75 | 0.125 | 0.111 | 0.098 | 0.082 | 0.017
56.25 | 0.151 | 0.130 | 0.120 | 0.105 | 0.012
75 | 0.209 | 0.189 | 0.178 | 0.163 | 0.011
$\mu$ (m2V-1s-1) | 12.5 | 0.84 | 3.71 | 6.68 | 1.67 | 4.77
25 | 0.46 | 0.78 | 0.48 | 0.85 | 3.34
43.75 | 0.28 | 0.30 | 0.34 | 0.40 | 1.97
56.25 | 0.22 | 0.26 | 0.28 | 0.32 | 2.78
75 | 0.16 | 0.18 | 0.19 | 0.20 | 3.04
Table 1: Cohesive energy, band gap, effective mass, and electron mobility for
various BN-doped graphene systems.

We next calculate the electron mobility for all the systems under study. The
electron mobility can be expressed as $\mu=\frac{e\tau}{m^{*}}$, where $\tau$
is the momentum relaxation time with a typical value of 1.9$\cdot 10^{-13}$ s
for Si doped graphene fedorov  and e is the charge of electron (1.6$\cdot
10^{-}{19}$ c). The estimated $\mu$ varies from 0.18 m2V-1s-1 to 6.68
m2V-1s-1, see Table I. We found the $\mu^{*}$ is amazingly high for the BN-
doped trilayer graphene which is 6 times higher than that of SiO2 gated
pristine graphene Ishigami  and 1.11 times higher than that of graphene on
$h$-BN substrate hone . We believe that our investigated systems should be
excellent candidate for device application.

In conclusion, we performed the comparative study of electronic properties of
BN doped monolayer graphene, bilayer, trilayer, multilayer, and superlattice
of pristine and BN-doped graphene. Our results prevail that the band gap of
graphene system can be tuned by BN-doping with different stacking pattern and
concentrations. Due to the breaking of the sublattice symmetry by presence of
B and N atoms a band gap of 0.02 eV to 2.43 eV is obtained for 12.5 % to 75%
BN-doping concentration. Our estimated electron effective mass for the systems
under study is found to be varies form 0.007 me to 0.209 me, which is smaller
than that of obtained for similar systems before. A giant mobility is
estimated as compared to experimentally observed value for graphene on $h$-BN
hone , would have great impact on devices. Finally, our predict is to open
finite and tunable band gap in BN-doped graphene with a giant mobility, which
is promising for future graphene based electronics ajayan ; ajayan1 ; Kim .

## References

  * (1) S. Y. Zhou, G.-H. Gweon, A. V. Fedorov, P. N. First, W. A. De Heer, D. -H. Lee, F. Guinea, A. H. Castro Neto, and A. Lanzara, Nat. Mater. 6, 770 (2007).
  * (2) T. P. Kaloni, Y. C. Cheng, and U. Schwingenschlögl, J. Mater. Chem. 22, 919 (2012).
  * (3) G. Giovannetti, P. A. Khomyakov, G. Brocks, P. J. Kelly, and J. van den Brink, Phys. Rev. B 76, 073103 (2007).
  * (4) Y. Xu, Y. Liu, H. Chen, X. Lin, S. Lin, B. Yu, and J. Luo, J. Mater. Chem. 22, 23821 (2012).
  * (5) L. Ci, L. Song, C. Jin, D. Jariwala, D. Wu, Y. Li, A. Srivastava, Z. F. Wang, K. Storr, L. Balicas, F. Liu, and P. M. Ajayan, Nat. Mater. 9, 430 (2010).
  * (6) C.-H. Park and S. G. Louie, Nano Lett. 10, 426 (2010).
  * (7) Y. W. Son, M. L. Cohen, and S. G. Louie, Phys. Rev. Lett. 97, 216803 (2010).
  * (8) G. Li, A. Luican, J. M. B. Lopes dos Santos, A. H. Castro Neto, A. Reina, J. Kong, and E. Y. Andrei, Nat. Phys. 6, 109 (2010).
  * (9) A. Bagri, C. Mattevi, M. Acik, Y. J. Chabal, M. Chhowalla, and V. B. Shenoy, Nat. Chem. 2, 581 (2010).
  * (10) J. Dai and J. Yuan, Phys. Rev. B 81, 165414 (2010).
  * (11) A. Nourbakhsh, M. Cantoro, T. Vosch, G. Pourtois, F. Clemente, M. H. van der Veen, J. Hofkens, M. M. Heyns, S. D. Gendt, and B. F. Sels, Nanotechnology 21, 435203 (2010).
  * (12) M. Pumera and C. H. A. Wong, Chem. Soc. Rev. 42, 5987 (2013).
  * (13) F. Withers, T. H. Bointon, M. Dubois, S. Russo, and M. F. Craciun, Nano Lett. 11, 3912 (2011).
  * (14) T. P. Kaloni, Y. C. Cheng, and U. Schwingenschlögl, EPL 100, 37003 (2013).
  * (15) T. P. Kaloni, M. Upadhyay-Kahaly, R. Faccio, and U. Schwingenschlögl, Carbon 64, 281 (2013).
  * (16) Z. Liu, L. Ma, G. Shi, W. Zhou, Y. Gong, S. Lei, X. Yang, J. Zhang, J. Yu, K. P. Hackenberg, A. Babakhani, J.-C. Idrobo, R. Vajtai, J. Lou, and P. M. Ajayan, Nat. Nanotechnol. 8, 119 (2013).
  * (17) K. Watanabe, T. Taniguchi, and H. Kanda, Nat. Mater. 3, 404 (2004).
  * (18) Y. Kubota, K. Watanabe, O. Tsuda, and T. Taniguchi, Science 317, 932 (2007).
  * (19) A. A. Balandin, S. Ghosh, W. Bao, I. Calizo, D. Teweldebrhan, F. Miao, and C. N. Lau, Nano. Lett. 8, 902 (2008).
  * (20) I. W. Frank, D. M. Tanenbaum, A. M. Van der Zande, and P. L. McEuen, J. Vac. Sci. Technol. 25, 2558 (2007).
  * (21) L. Song, L. Ci, H. Lu, P. B. Sorokin, C. Jin, J. Ni, A.G. Kvashnin, D. G. Kvashnin, J. Lou, B. I. Yakobson, and P. M. Ajayan, Nano. Lett. 10, 3209 (2010).
  * (22) J. G. Kho, K. T. Moon, J. H. Kim, and D. P. J. Kim, Am. Ceram. Soc. 83, 2681 (2000).
  * (23) Y. Chen, J. Zou, S. J. Campbell, and G. Le Caer, Appl. Phys. Lett. 84, 2430 (2004).
  * (24) J. Han, L. L. Zhang, S. Lee, J. Oh, K.-S. Lee, J. R. Potts, J. Ji, X. Zhao, R. S. Ruoff, and S. Park, ACS Nano 7, 19 (2013).
  * (25) T. Sharifi, E. Gracia-Espino, H. R. Barzegar, X. Jia, F. Nitze, G. Hu, P. Nordblad, C.-W. Tai, and T. Wágberg, Nat. Commun. 4, 2319 (2013).
  * (26) Y. Xue, D. Yu, L. Dai, R. Wang, D. Li, A. Roy, F. Lu, H. Chen, Y. Liu, and J. Qu, Phys. Chem. Chem. Phys. 15, 12220 (2013).
  * (27) T. P. Kaloni and S. Mukherjee, Mod. Phys. Lett. B 25, 1855 (2011).
  * (28) O. Leenaerts, H. Sahin, B. Partoens, and F. M. Peeters, Phys. Rev. B 88, 035434 (2013).
  * (29) T. P. Kaloni and S. Mukhergie, J. Nanopart. Res. 14, 1059 (2012)
  * (30) S. M. Kim, A. Hsu, P. T. Araujo, Y.-H. Lee, T. Palacios, M. Dresselhaus, J.-C. Idrobo, K. K. Kim, and J. Kong, Nano Lett. 13, 933 (2013).
  * (31) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1986).
  * (32) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiaro tti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquare llo, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, J. Phys. Condens. Matt. 21, 395502 (2009).
  * (33) S. Grimme, J. Comput. Chem. 27, 1787 (2006).
  * (34) N. N. Greenwood and A. Earnshaw, Chemistry of the Elements (Pergamon, Oxford, 1984).
  * (35) M. Topsakal, E. Aktürk, and S. Ciraci, Phys. Rev. B 79, 115442 (2009).
  * (36) C. Chang, S. Kataria, C.-C. Kuo, A. Ganguly, B.-Y. Wang, J.-Y. Hwang, K.-J. Huang, W.-H. Yang, S.-B. Wang, C.-H. Chuang, M. Chen, C.-I. Huang, W.-F. Pong, K.-J. Song, S.-J. Chang, J.-H. Guo, Y. Tai, M. Tsujimoto, S. Isoda, C.-W. Chen, L.-C. Chen, and K.-H. Chen, ACS Nano 7, 1333 (2013).
  * (37) L. Liu and Z. Shen, Appl. Phys. Lett. 95, 252104 (2009).
  * (38) M. Freitag, Nat. Phys. 7, 596 (2011).
  * (39) Y. Zhang, T.-T. Tang, C. Girit, Z. Hao, M. C. Martin, A. Zettl, M. F. Crommie, Y. R. Shen, and F. Wang, Nature 459, 820 (2009).
  * (40) B. Xu, Y. H. Lu, Y. P. Feng, and J. Y. Lin, J. Appl. Phys. 108, 073711 (2010).
  * (41) S. Casolo, R. Martinazzo, and G. F. Tantardini, J. Phys. Chem. C 115, 3250 (2011).
  * (42) Y. Fan, M. Zhaoa, Z. Wang, X. Zhang, and H Zhang, Appl. Phys. Lett. 98, 083103 (2011).
  * (43) R. Magno, E. R. Glaser, B. P. Tinkham, J. G. Champlain, J. B. Boos, M. G. Ancona, and P. M. Campbell, J. Vac. Sci. Technol. B 24, 1622 (2006).
  * (44) D. V. Fedorov, M. Gradhand, S. Ostanin, I. V. Maznichenko, A. Ernst, J. Fabian, and I. Mertig, Phys. Rev. Lett. 110, 156602 (2013).
  * (45) J.-H. Chen, C. Jang, M. Ishigami, S. Xiao, W. G. Cullen, E. D. Williamsa, M. S. Fuhrera, Solid State Commun. 149, 1080 (2009).
  * (46) C. R. Dean, A. F. Young, I. Meric, C. Lee, L. Wang, S. Sorgenfrei, K. Watanabe, T. Taniguchi, P. Kim, K. L. Shepard, and J. Hone, Nat. Nanotechnol. 5, 722 (2010).