# GenMarkov: Modeling Generalized Multivariate Markov Chains in R

Carolina Vasconcelos, Bruno Damásio
NOVA Information Management School (NOVA IMS)
Universidade Nova de Lisboa
Campus de Campolide, 1070-312 Lisboa, Portugal
{cvasconcelos<EMAIL_ADDRESS>

###### Abstract

This article proposes a new generalization of the Multivariate Markov Chains
(MMC) model. The future values of a Markov chain commonly depend on only the
past values of the chain in an autoregressive fashion. The generalization
proposed in this work also considers exogenous variables that can be
deterministic or stochastic. Furthermore, the effects of the MMC’s past values
and the effects of pre-determined or exogenous covariates are considered in
our model by considering a non-homogeneous Markov chain. The Monte Carlo
simulation study findings showed that our model consistently detected a non-
homogeneous Markov chain. Besides, an empirical illustration demonstrated the
relevance of this new model by estimating probability transition matrices over
the space state of the exogenous variable. An additional and practical
contribution of this work is the development of a novel R package with this
generalization.

_K_ eywords Multivariate Markov chains and Mixture transition distribution
Model and High order Markov chains and Multivariate Markov chains with
exogenous variables and R

##  1 Introduction

Multivariate Markov chains (MMC) have a wide range of applications, in various
fields. Hence, several studies and generalizations of the MMC models have been
made. However, the avail- ability of packages that allow the estimation and
application of these models are scarce, and most of these methods use
algorithms and software that are not broadly available or can only be applied
in particular situations.

In the last few years, R software has been gaining importance in the field of
statistical computing. This phenomenon might be because it is free and open-
source software, which compiles and runs on a wide variety of operating
systems.

Specifically, in R software, there are some available packages related to
Markov chains (MC) and MMC. For example, the march package (Maitre, Emery,
2020; Berchtold et al., 2020) allows the computation of various Markovian
models for categorical data, including homogeneous Markov chains of any order,
MTD models, Hidden Markov models, and Double Chain Markov Models. Ogier Maitre
developed this package with contributions from Andre Berchtold, Kevin Emery,
Oliver Buschor, and Andre Berchtold maintains it. All the models computed by
this package are for univariate categorical data. The markovchain package
(Spedicato, 2017) contains functions and methods to create and manage
discrete-time Markov chains. In addition, it includes functions to perform
statistical and probabilistic analysis (analysis of their structural
proprieties). Finally, the DTMCPack package (Nicholson, 2013) contains a
series of functions that aid in both simulating and determining the properties
of finite, discrete-time, discrete-state Markov chains. There are two main
functions: DTMC and MultDTMC, which produce $n$ iterations of a Markov
Chain(s) based on transition probabilities and an initial distribution given
by the user, for the univariate and multivariate case, respectively. This last
package is the only one available in R for MMC.

The main goal of this paper is not only the development of the Generalized
Multivariate Markov chain (GMMC) models and addressing statistical inference
in MMC models, but also the implementation of these methods in an R package.
The R package includes three functions: multimtd, multimtd_probit and mmcx.
The first two functions estimate the MTD model for multivariate categorical
data, with Chings’s specification (Ching et al., 2002) and with the Probit
specification (Nicolau, 2014), respectively. The last function allows the
estimation of our proposed model, the Generalized Multivariate Markov Chain
(GMMC) model. The R package, GenMarkov, with these three functions is
available in the Comprehensive R Archive Network (CRAN) at
https://CRAN.R-project.org/package=GenMarkov.

##  2 Multivariate Markov Chains

Markov chains can be appropriate for representing dependencies between
successive observations of a random variable. However, when the order of the
chain or the number of possible values increases, Markov chains have lack
parsimony. In this context, Jacobs, Lewis (1978); Pegram (1980); Logan (1981)
proposed several models for HOMC. Notwithstanding these developments, the
Mixture Transition Distribution model (Raftery, 1985) proved to be more
suitable to model HOMC, which overshadowed the previously proposed models.

Several relevant extensions of the MTD model emerged: the Multimatrix MTD
(Berchtold, 1995, 1996), which allowed modeling the MTD by using a different
$m\times m$ transition matrix for each lag, the Infinite-Lag MTD model that
assumes an infinite lag order ($l=\infty$), which was first considered by
Mehran (1989) and later developed by Le et al. (1996) in a more general
context. Finally, the MTD with General State Spaces allowed modeling more
general processes with an arbitrary space state (Martin, Raftery, 1987; Adke,
Deshmukh, 1988; Wong, Li, 2001).

Although the MTD model presents a more parsimonious approach to model Markov
chains with order higher than one, it has weaknesses. Namely, when considering
more than one data sequence, one represents the MMC as a HOMC, by expanding
the state-space. This approach could result in a more complex probability
transition matrix. Consequently, this can make the estimation unfeasible as
the order, states, and the number of data sequences increase. Additionally,
the model assumes the same transition matrix for each lag.

In this setting, Ching et al. (2002) determined an alternative to handle the
unfeasibility of the conventional multivariate Markov chain (MMC) by proposing
a model with fewer parameters. The model developed is essentially the same as
the MTD. However, it considers a different $m\times m$ transition matrix for
each lag and considers more than one data sequence.

In the proposed multivariate Markov chain model, Ching et al. (2002) assume
the following relationship:

Let $x_{t}^{(j)}$ be the state vector of the $j$th sequence at time $t$. If
the $j$th sequence is in state $l$ at time $t$ then

$x_{t}^{(j)}=(0,\dots,0,\underbrace{1}_{\text{jth entry}},0,\dots,0)^{t}$

The model is given by

$x_{t+1}^{(j)}=\sum_{k=1}^{s}\lambda_{jk}P^{(jk)}x_{t}^{(k)},\text{for
}j=1,2,\dots,s$ (1)

where $\lambda_{jk}\geq 0$ for $1\leq j,k\leq s$ and
$\sum_{k=1}^{s}\lambda_{jk}$ for $j=1,2,\dots,s$.

The state probability distribution of the $k$th sequence at time $(t+1)$
depends on the weighted average of $P^{(jk)}x_{t}^{(k)}$ . Here $P^{(jk)}$ is
a transition probability matrix from the states in the $k$th sequence to the
states in the $j$th sequence and $x_{t}^{(k)}$ is the state probability
distribution of the $k$th sequences at time $t$. In matrix form:

$\underline{x}_{t+1}^{(j)}\equiv\left[\begin{array}[]{c}x_{t+1}^{(1)}\\\
\vdots\\\
x_{t+1}^{(s)}\end{array}\right]=\left[\begin{array}[]{ccc}\lambda_{11}P^{(11)}&\dots&\lambda_{1s}P^{(1s)}\\\
\vdots&\ddots&\vdots\\\
\lambda_{s1}P^{(s1)}&\dots&\lambda_{ss}P^{(ss)}\end{array}\right]\left[\begin{array}[]{c}x_{t}^{(1)}\\\
\vdots\\\ x_{t}^{(s)}\end{array}\right]\equiv Q\underline{x}_{t}$ (2)

where $Q$ is an $ms\times ms$ block matrix ($s\times s$ blocks of $m\times m$
matrices) and $x_{t}$ is a stacked $ms$ column vector ($s$ vectors, each one
with $m$ rows).

The matrices $P^{(jk)}$ can be estimated for each data sequence by counting
the transition frequency from the states in the $k$th sequence those in the
$j$th sequence, obtaining the transition frequency matrix for the data
sequence. After normalization, the estimates of the transition probability
matrices, i.e., $\hat{P}^{(jk)}$, are obtained.

Regarding the $\lambda_{jk}$ coefficients, the estimation method proposed by
Ching et al. (2002) involves the following optimization problem:

$\displaystyle\begin{split}min_{\lambda}max_{i}|[\sum_{k=1}^{m}\lambda_{jk}\hat{P}^{(jk)}\hat{\boldsymbol{x}}^{(k)}-\hat{\boldsymbol{x}}^{(j)}]|\\\
\text{s.t. }\sum_{k=1}^{s}\lambda_{jk}\text{ and }\lambda_{jk}\geq
0\end{split}$ (3)

Besides this, different models have been proposed for multiple categorical
data sequences. Kijima et al. (2002) proposed a parsimonious MMC model to
simulate correlated credit risks. Siu et al. (2005) proposed an easy to
implement model; however, its applicability was limited by the number of
parameters involved. Ching et al. (2008) proposed a simplified model based on
an assumption proposed in Zhang et al. (2006). Zhu, Ching (2010) proposed a
method of estimation based on minimizing the prediction error with equality
and inequality restrictions and Nicolau, Riedlinger (2014) proposed a new
approach to estimate MMC which avoids imposing restrictions on the parameters,
based on non-linear least squares estimation, facilitating the model
estimation and the statistical inference. (Berchtold, 2003) proposed a MTD
model for heteroscedastic time series. Lastly, Wang et al. (2014) proposed a
new multivariate Markov chain model to reduce the number of parameters. Thus,
generally, the models used in the published papers were developed by Ching et
al. (2002) or were a consequent generalization of them and addressed the MMC
as an end in itself.

In Damásio (2013); Damásio, Nicolau (2014), a different and innovative concept
was proposed: the usage of MMC as regressors in a certain model. Hence, given
that the MMC Granger causes a specific dependent variable, and taking
advantage of the information about the past state interactions between the MMC
categories, it was possible to forecast the current dependent variable more
accurately.

Other relevant contributions are related to the optimization algorithm, as in
Lèbre, Bourguignon (2008) and Chen, Lio (2009), and to empirical applications
(Ching et al., 2003; Ching, Ng, 2006; Damásio et al., 2018; Damásio, Mendonça,
2019, 2020). Also, Damásio, Nicolau (2020) proposed a new methodology for
detecting and testing the presence multiple structural breaks in a Markov
chain occurring at unknown dates.

In the vast majority of MMC models’ studies, a positive correlation between
the different data sequences is assumed due to the restrictions imposed. This
aspect means it is always considered that at moment $t$, an increase in a
state probability for a data sequence has an increasing impact on another data
sequence, for time $t+1$. Thereupon, if one has a negative correlation between
series, the parameter estimates are forced to be zero. The solution to this
problem is very straightforward; one can relax the assumptions and not assume
the constraints. However, that means the results produced by the model will no
longer be probabilities. Raftery, Tavaré (1994) presented an alternative, by
dropping the positivity condition and imposing another set of restrictions.
Ching et al. (2008) also tackled this issue and proposed a method where one
splits the $Q$ matrix into the sum of two other matrices and one represents
the positive correlations and another the negative correlations. Also, in
Nicolau (2014), a specification completely free from constraints, inspired by
the MTD model, was proposed, facilitating the estimation procedure and, at the
same time, providing a more accurate specification for
$P_{j}(i_{0}|i_{1},\dots,i_{s})$. The model was:

$P_{j}(i_{0}|i_{1},\dots,i_{s})=P_{j}^{\Phi}(i_{0}|i_{1},\dots,i_{s}):=\\\
\frac{\Phi(\eta_{j0}+\eta_{j1}P(i_{0}|i_{1})+\dots+\eta_{js}P(i_{0}|i_{s}))}{\sum_{k=1}^{m}\Phi(\eta_{j0}+\eta_{j1}P(k|i_{1})+\dots+\eta_{js}P(k|i_{s}))}$
(4)

where $n_{ji}\in\mathbb{R}(j=1,\dots,s;i=1,\dots,m)$ and $\Phi$ is the
(cumulative) standard normal distribution function. This specification is
denoted as and MTD-Probit model. The log-likelihood is given by:

$LL=\sum_{i_{1},i_{2},\dots,i_{i_{s}},i_{0}}n_{i_{1},i_{2},\dots,i_{i_{s}},i_{0}}log(P_{j}^{\Phi}(i_{0}|i_{1},\dots,i_{s}))$
(5)

and the maximum likelihood estimator is defined, as usual, as
$\hat{\eta}=\text{arg max}_{n_{j1},\dots,n_{js}}LL$. The parameters
$P_{jk}(i_{0}|i_{1})$, $k$ =$1,\dots,s$ can be estimated in advance, through
the consistent and unbiased estimators proposed by Ching et al. (2002):

$\widehat{P_{jk}(i_{0}|i_{1})}=\frac{n_{i_{1}i_{0}}}{\sum_{i_{0}=1}^{n}n_{i_{1}i_{0}}}$
(6)

This specification can be superior to the MTD because the estimation procedure
is easier, and the standard numerical optimization routines can be easily
applied in the absence of constraints. However, similarly to the standard MTD,
the likelihood is not a strictly concave function on the entire parameter
state-space, thus the choice of starting values is still important.
Additionally, the model describes a broader range of possible dependencies
since the parameters are not constrained. Moreover, this proposed model is
more accurate than the MTD model. For more details on this, see Nicolau
(2014).

Overall, the published work on MMC models was mostly based on improving the
estimation methods and/or making the model more parsimonious. In Damásio
(2013); Damásio, Nicolau (2014), a different approach was used, and the work
developed focused on the usage of MMC as regressors in a certain model.
Notably, it showed that an MMC can improve the forecast of a dependent
variable. In a way, it demonstrated that an MMC can be an end in itself, but
it can be an instrument to reach an end or a purpose. In this work, the
opposite will be developed: instead of considering an MMC as regressors, a
model in which a vector with pre-determined exogenous variables is part of
$\mathcal{F}_{t-1}$ is proposed.

##  3 Covariates in Markov Chain Models

Regarding the inclusion of covariates in Markov chains models, Regier (1968)
proposed a two-state Markov chain model, where the transition matrix
probabilities were a function of a parameter, $q$, that described the tendency
of the subject to move from state to state.Kalbfleisch, Lawless (1985)
proposed a panel data analysis method under a continuous-time Markov model
that could be generalized to handle covariate analysis and the fitting of
certain non-homogeneous models. This work overcame the limitations of
Bartholomew (1968), Spilerman, Singer (1976) and Wasserman (1980)
methodologies, by developing a new algorithm that provided a very efficient
way of obtaining maximum likelihood estimates. Also, Muenz, Rubinstein (1985)
developed a Markov model for covariates dependence of binary sequences, where
the transitions probabilities were estimated through two logistic regressions
that depended on a set of covariates. Essentially, Muenz, Rubinstein (1985)
modeled a non-homogeneous Markov chain through logistic regression,
considering only two states. Islam et al. (2004) developed an extension of
this model considering three states, and Islam, Chowdhury (2006) generalized
this approach for HOMC. Additionally, Azzalini (1994) proposed a model to
study the influence of time-dependent covariates on the marginal distribution
of a binary response in serially correlated binary data, where Markov chains
are expressed in terms of transitional probabilities.

More recently, Bolano (2020) proposed an MTD-based approach to handle
categorical covariates, that considers each covariate separately and combines
the effects of the lags of the MTD and the covariates employing a mixture
model. Specifically, the model is given by:

$P(X_{t}=k\mid
X_{t-1}=i,C_{1}=c_{1},\dots,C_{l}=c_{l})\approx\theta_{0}a_{ik}+\sum_{h=1}^{l}\theta_{h}d_{c_{h}k}$
(7)

where $a_{ik}$ is the transition probability from state $i$ to state $k$, as
in a conventional Markov chains and $d_{c_{h}k}$ is the probability of
observing the states $k$ given the modality $c_{h}$ of the covariate $h$.
Lastly, $\theta_{0},\dots,\theta_{l}$ are the weights of the explanatory
elements of the model.

According to the literature presented, several researchers have proposed
methodologies or generalizations to include covariates in Markov chain models.
Primarily for social sciences and health applications, where the transition
probabilities were generally modeled through logistic regression. However,
there has been an increased focus on categorical covariates, opposing
continuous covariates and a lack of approaches to multivariate Markov chain
models. Thus, with this work, we aim to tackle this research gap.

##  4 Multivariate Markov Chains with covariates

###  4.1 Theoretical model

In this work, a new generalization of Ching et al. (2002) MMC model is
presented: the GMMC model, that is, we will consider exogeneous or pre-
determined covariates in the $\sigma$ \- algebra generated by the available
information until $t-1$ ($\mathcal{F}_{t-1}$). These variables can be
deterministic or stochastic and do not necessarily need to be reported at time
$t$. Broadly, the model is given by:

$P(S_{jt}=k|\mathcal{F}_{t-1})=P(S_{jt}=k|S_{1t-1}=i_{1},S_{2t-1}=i_{2},\dots,S_{st-1}=i_{s},\boldsymbol{x}_{t})$
(8)

We can specify this model as proposed by Ching et al. (2002) with Raftery’s
notation:

$P(S_{jt}=i_{0}|S_{1t-1}=i_{1},\dots,S_{st-1}=i_{s},\boldsymbol{x}_{t})\equiv\\\
\lambda_{j1}P(S_{jt}=i_{0}|S_{1t-1}=i_{1},\boldsymbol{x}_{t})+\dots+\lambda_{js}P(S_{jt}=i_{0}|S_{st-1}=i_{s},\boldsymbol{x}_{t})$
(9)

subject to the usual constraints.

###  4.2 Estimation and Inference

This proposed model is estimated through MLE, similar to the standard MTD
model. The log-likelihood is given by:

$LL=\sum_{t=1}^{n}logP(S_{jt}=i_{0}|S_{1t-1}=i_{1},\dots,S_{st-1}=i_{s},\boldsymbol{x}_{t})$
(10)

Additionally, the probabilities can be estimated through an multinomial logit
model. The following Proposition holds the consistency of the MLE estimator:

######  Proposition (2).

Let $\\{\boldsymbol{w}_{t}\\}$ be ergodic stationary random variable, with
likelihood function $f(\boldsymbol{w}_{t}|\theta_{0})$. Let $\hat{\theta}$ be
the MLE estimator. Suppose that:

  1. (i)

$E[\mbox{log }f(\boldsymbol{w}_{t}\mid\theta_{0})]$ is uniquely maximized on
$\Theta$ at $\theta_{0}\in\Theta$,

  2. (ii)

$\theta_{0}\in\Theta$, which is compact,

  3. (iii)

$\mbox{log }f(\boldsymbol{w}_{t}\mid\theta_{0})$ is continuous at each
$\theta\in\Theta$ with probability one,

  4. (iv)

$E[sup_{\theta\in\Theta}\mid\mbox{log }f(w_{t}\mid\theta)\mid]<\infty$

Then $\hat{\theta}\xrightarrow{p}\theta_{0}$

Condition (i) is verified according to Lemma 2.2 of Newey, Mcfadden (1994).
Condition (ii) is verified and guaranteed by the restrictions imposed in the
model parameters. Knowing that
$P(S_{jt}=i_{0}|S_{1t-1}=i_{1},\dots,S_{st-1}=i_{s},\boldsymbol{x}_{t})$ is
linear combination of a set of $n$ probabilities and since the logarithm
function is a continuous function, condition (iii) is verified. Finally,
condition (iv) is verified according to Lemma 2.4 of Newey, Mcfadden (1994).

Regarding inference, MLE will be asymptotically normal if it is consistent and
the following Proposition verifies:

######  Proposition (3).

Let $\\{\boldsymbol{w}_{t}\\}$ be ergodic stationary random variable and let
$s(\boldsymbol{w}_{t};\theta)$ and $H(\boldsymbol{w}_{t};\theta)$ be the first
and second partial derivatives of the $\mbox{log
}f(\boldsymbol{w}_{t}\mid\theta)$, respectively. Suppose the estimator
$\hat{\theta}$ is consistent and suppose, further, that

  1. (i)

$\theta_{0}$ is in the interior of $\Theta$,

  2. (ii)

$\mbox{log }f(\boldsymbol{w}_{t}\mid\theta_{0})$ is twice continuously
differentiable in $\theta$ for any $\boldsymbol{w}_{t}$,

  3. (iii)

$\frac{1}{\sqrt{n}}\sum^{n}_{t=1}s(\boldsymbol{w}_{t};\theta_{0})\underset{d}{\rightarrow}N(0,\Sigma)$,
where $\Sigma$ is positive definite,

  4. (iv)

For some neighborhood of $\mathcal{N}$ of $\theta_{0}$,

$E[\underset{\theta\in\mathcal{N}}{sup}||H(\boldsymbol{w}_{t};\theta)||]<\infty$

so that for any consistent estimator $\tilde{\theta}$,
$\frac{1}{n}\sum_{t=1}^{n}H(\boldsymbol{w}_{t};\tilde{\theta})\underset{p}{\rightarrow}E[H(\boldsymbol{w}_{t};\theta)]$

  5. (v)

$E[H(\boldsymbol{w}_{t},\theta_{0})]$ is nonsingular.

Then $\hat{\theta}$ is asymptotically normal with

$Avar(\hat{\theta})=\\{E[H(\boldsymbol{w}_{t};\theta_{0})]\\}^{-1}\Sigma\\{E[H(\boldsymbol{w}_{t};\theta_{0})]\\}^{-1}$

Condition (i) is verified as in condition (ii) of Proposition (2). Condition
(ii) is also verified as in condition (iii) of Proposition (2), since the
logarithm function is twice continuously differentiable. Condition (iii) is
verified according to the Ergodic Stationary Martingale Differences CLT
(Billingsley, 1961). In this case, $\Sigma$ =
$E[s(\boldsymbol{w}_{t};\theta_{0})s(\boldsymbol{w}_{t};\theta_{0})^{\prime}]$
= $-E[H(\boldsymbol{w}_{t};\theta_{0})]$, which implies that
$Avar(\hat{\theta})=-\\{E[H(\boldsymbol{w}_{t};\theta_{0})]\\}^{-1}$.
Condition (iv) is verified according to Lemma 2.4 of Newey, Mcfadden (1994).
Considering only one equation, let $\boldsymbol{q}_{t}$ be a $t\times s$
matrix of the probabilities $P(S_{1t}\mid S_{1t},x_{t}),\dots,P(S_{1t}\mid
S_{st},x_{t})$ and $\boldsymbol{\lambda}$ a row-vector of
$\lambda_{11},\dots,\lambda_{1s}$, the hessian matrix is given by
$E[\boldsymbol{q}_{t}^{\prime}\boldsymbol{q}_{t}[(\boldsymbol{\lambda}\boldsymbol{q}_{t}^{\prime})(\boldsymbol{\lambda}\boldsymbol{q}_{t}^{\prime})^{\prime}]^{-1}]$.
Condition (v) is verified if
$E[\boldsymbol{q}_{t}^{\prime}\boldsymbol{q}_{t}]$ is nonsingular.

###  4.3 Monte Carlo Simulation Study

A Monte Carlo simulation study was designed to evaluate the dimension and
power of the test parameters of the proposed model. The R statistical
environment was used for all computations. This simulation study was comprised
of two parts.

####  4.3.1 Part I: Detect a non-homogeneous Markov chain

First, we considered two sequences with two and three states. The main goal
was to assess if the model detected the presence of a non-homogeneous Markov
chain correctly and if the estimate of the parameter would correspond to the
expected. So, given two sequences, one generated through a non-homogeneous
Markov chain and the other generated through a homogeneous Markov chain, it
would be expected that the parameter associated with the transition
probabilities of the first sequence would be one and the parameter associated
with the transition probabilities of the second sequence would be zero. With
this in mind, the transitions probabilities of the first sequence were
estimated through a logistic regression, where parameters of this regression
were randomly generated in R, and the second sequence was generated through a
first-order Markov chain. Hence, for both states cases considered, it was
expected that the estimated regression would be:

$P(S_{1t}=i_{0}|S_{1t-1}=i_{1},S_{2t-1}=i_{2},\boldsymbol{x}_{t-1})=\\\
1\times P(S_{1t}=i_{0}|S_{1t-1}=i_{1},\boldsymbol{x}_{t-1})+0\times
P(S_{1t}=i_{0}|S_{2t-1}=i_{2},\boldsymbol{x}_{t-1})$ (11)

To assess the test power and dimension, we used the Wald test with the
following hypothesis:

| Hypothesis | Test
---|---|---
Power | $H_{0}:\lambda_{11}=0$ | $\frac{\hat{\lambda}_{11}^{2}}{se(\hat{\lambda}_{11})^{2}}\sim\chi^{2}_{(1)}$
$H_{0}:\lambda_{12}=1$ | $\frac{(\hat{\lambda}_{12}-1)^{2}}{se(\hat{\lambda}_{12})^{2}}\sim\chi^{2}_{(1)}$
Dimension | $H_{0}:\lambda_{11}=1$ | $\frac{(\hat{\lambda}_{11}-1)^{2}}{se(\hat{\lambda}_{11})^{2}}\sim\chi^{2}_{(1)}$
$H_{0}:\lambda_{12}=0$ | $\frac{\hat{\lambda}_{12}^{2}}{se(\hat{\lambda}_{12})^{2}}\sim\chi^{2}_{(1)}$

The simulation procedure was performed as follows:

  1. 1.

Generate the values of the coefficients for the probability transition matrix
of series $S_{1t}$ randomly;

  2. 2.

Generate the probability transition matrix of series $S_{2t}$ randomly;

  3. 3.

Set the initial value of $S_{2t}$ to 1 and simulate the following from the
defined probability transition matrix;

  4. 4.

In each iteration (of 1000 repetitions),

     * •

Generate $X_{t}\sim N(2,25)$;

     * •

Generate the time-varying probabilities of series $S_{1t}$ through the values
of the fixed coefficients and the lagged variable $x_{t}$;

     * •

Set the initial values of the series $S_{1t}$ as 1;

     * •

For each period $t$, simulate the next state of $S_{1t}$ from the
probabilities simulated for that moment;

     * •

Estimate the model through the function mmcx;

     * •

Calculate the Wald test and add to the counter if it is rejected.

(a) Test Dimension

(b) Test Power

Figure 1: Simulation study results: Two-states

Considering two states, the test dimension was at 5.7% with a sample size of
100 observations, sightly increased with 500 observations, and returned to the
expected values in 1000 and 5000 observations. For a sample size of 100, 500,
and 1000 observations, we have low test power. So, when considering two
states, the sample must have at least 5000 observations, or, if that is not
possible, consider a higher significance level when testing for individual
significance.

(a) Test Dimension

(b) Test Power

Figure 2: Simulation study results: Three-states

Considering three states, the test dimension was 9.7% for a sample size of 100
observations, 0.2% for a sample size of 500 observations, and 0.3% for a
sample size of 1000. Regarding the test power, we see similar behavior, for a
sample of 100 observations, the test power was 90.5%, and from a sample of 500
observations, we reach a test power of 100%. Thus, when considering three
states, one may consider a sample of 500 observations without compromising the
test power and dimension.

####  4.3.2 Part II: Detect parameters assigned values

Secondly, we performed a simulation study where we considered two non-
homogeneous Markov chains with two states. Here, the main goal was to assess
if the model correctly detected the parameters assigned. So, in this case, we
started by generating the terms of the model proposed. These terms were
estimated through logistic regression, and the parameters of this regression
were randomly generated in R. Similarly to Part I, we considered a Wald test
to assess the power and dimension of the test. The simulation procedure was
performed as follows:

  1. 1.

Generate the values of the coefficients to calculate the probability
transition matrices randomly;

  2. 2.

In each iteration (of 1000 repetitions),

     * •

Generate $X_{t}\sim N(2,25)$;

     * •

Generate the probabilities $P(S_{jt}|S_{st-1},x_{t-1})$, with $j=1,2$ and
$s=1,2$.

     * •

Set the initial values of the series $S_{1t}$ and $S_{2t}$ as 1;

     * •

For each period $t$, calculate the probabilities
$P(S_{1t}|S_{1t-1},S_{2t-1},x_{t-1})$ and $P(S_{2t}$
$|S_{1t-1},S_{2t-1},x_{t-1})$ through the assigned values of the $\lambda$’s.
Considering the calculated probabilities, simulate the next state for each
series, $S_{1t}$ and $S_{2t}$.

     * •

Estimate the model through the function mmcx;

     * •

Calculate the Wald test and add to the counter if it is rejected.

The probabilities $P(S_{1t}|S_{1t-1},x_{t-1})$ and
$P(S_{1t}|S_{2t-1},x_{t-1})$ presented some differences regarding its values’
distributions. Specifically, $P(S_{1t}|S_{1t-1},x_{t-1})$ had more extreme
probabilities values, with the minimum value being close to 0 and the maximum
value being close to 1. And, the probabilities $P(S_{1t}|S_{2t-1},x_{t-1})$
had more moderate values, with the minimum value being, on average, 0.3 and
the maximum value, 0.7. When the probabilities have values close to 1, one
says that the states/regimes are persistent.

We calculated the power and dimension of test for each value of $\lambda$ when
the estimated probabilities are moderate and when they are extreme. Hence,
considering equation 1:

$P(S_{1t}=i_{0}|S_{1t-1}=i_{1},\dots,S_{2t-1}=i_{2},\boldsymbol{x}_{t-1})=\\\
\lambda_{11}P(S_{1t}=i_{0}|S_{1t-1}=i_{1},\boldsymbol{x}_{t-1})+\lambda_{2s}P(S_{1t}=i_{0}|S_{2t-1}=i_{s},\boldsymbol{x}_{t-1})$
(12)

The parameter $\lambda_{11}$ will be associated with more extreme
probabilities and $\lambda_{12}$ will be associated with more moderate
probabilities.

(a) Test Dimension

(b) Test Power

Figure 3: Simulation study results: Case 1 - Persistent states on low values
of the parameters

(a) Test Dimension

(b) Test Power

Figure 4: Simulation study results: Case 2 - Persistent states on high values
of the parameters

When the states are persistent and the parameter’s value is low (i.e., 0.2 and
0.4), we have low test power. By increasing this value, the power of test
increases as well. When the states are not persistent, we do not have a clear
pattern regarding the power of test, for a value of the parameter of 0.2, the
power of test is still low (although not as low as the first scenario),
increases when we have a value of 0.4, decreases when the value is 0.6 and
increases again when the value is 0.8. Overall, the estimated standard errors
seem high, leading to low test power. Regarding the test dimension, when we
have a higher weight associated with the non-persistent states, the test
dimension converges to 0. However, when this weight is associated with the
persistent states, the test dimension increases with the sample size, reaching
a value of 10% in some cases. Hence, one must use a 10% significance level to
perform statistical inference on the parameters in this situation.

###  4.4 Software Implementation

Regarding the software implementation for each function, for the multimtd
function the estimation method was presented in Berchtold (2001) applied to
the multivariate case. For multimtd_probit, a package for numerical
maximization of the log-likelihood, maxLik (Henningsen, Toomet, 2011), was
used. This package performs Maximum Likelihood estimation through different
optimization methods that the user can choose. The optimization methods
available are Newton-Raphson, Broyden - Fletcher - Goldfarb - Shanno, BFGS al-
algorithm, Berndt - Hall - Hall - Hausman, Simulated ANNealing, Conjugate
Gradients, and Nelder-Mead. Finally, for the mmcx function, a different
approach was used. Unlike the MTD- Probit, the model proposed has equality and
inequality restrictions in the parameters. The maxLik package only allows one
type of restriction for each Maximum Likelihood estimation, so it was not
possible to use this package to estimate the proposed model with exogenous
variables. Hence, the algorithm used was the Augmented Lagrangian method,
available in the alabama package through the function auglag. This estimation
method for the proposed model is not very common, however, it has been applied
to Markov chain models (Rajarshi, 2013). The GMMC model’s probabilities were
estimated through a Multinomial Logit using rmultinom of the neet package
(Venables, Ripley, 2002).

Additionally, the hessian matrices were also computed, which allowed
performing statistical inference. The maxLik and auglag compute the Hessian
matrices with the estimates. For the function multimtd, since the optimization
procedure of Berchtold (2001) was used, the hessian was computed through the
second partial derivatives.

The function multi.mtd requires the following elements:

  * •

y, a matrix of the categorical data sequences.

  * •

deltaStop, the delta below which the optimization phases of the parameters
stop.

  * •

is_constrained, flag indicating whether the function will consider the usual
set of constraints (usual set: TRUE, new set of constraints: FALSE).

  * •

delta, the amount of change to increase/decrease in the parameters for each
iteration of the optimization algorithm.

The last three arguments concern the optimization procedure. For more details
see Berchtold (2001).Considering two vectors of two categorical data
sequences, s1 and s2, to estimate the model and obtain the results:

    R> multi.mtd(y=cbind(s1,s2), deltaStop=0.0001, is_constrained=TRUE, delta=0.1)

The function multi.mtd_probit requires the following arguments:

  * •

y, a matrix of the categorical data sequences.

  * •

initial, a vector of the initial values of the parameters.

  * •

nummethod, the numerical maximization method, currently either "NR" (for
Newton-Raphson), "BFGS" (for Broyden-Fletcher-Goldfarb-Shanno), "BFGSR" (for
the BFGS algorithm implemented in R), "BHHH" (for Berndt-Hall-Hall-Hausman),
"SANN" (for Simulated ANNealing), "CG" (for Conjugate Gradients), or "NM" (for
Nelder-Mead). Lower-case letters (such as "nr" for Newton-Raphson) are
allowed. The default method is "BFGS". For more details see maxLik package.

Considering two vectors of two categorical data sequences, s1 and s2 again, to
estimate the model an obtain the results with BFGS maximization method:

    R> multi.mtd_probit(y = cbind(s1,s2), initial=c(1,1,1), nummethod=’bfgs’)

Finally, the function mmcx requires the following elements:

  * •

y, a matrix of categorical data sequences.

  * •

x, a matrix of covariates (exogeneous variables).

  * •

initial, a vector of the initial values of the parameters.

Considering two vectors of two categorical data sequences, s1 and s2, and a
vector of an exogeneous variables, x, to estimate the model and obtain the
results:

    R> mmcx(y = cbind(s1,s2), x = cbind(x), initial=c(1,1))

These functions return a list with the parameter estimates, standard errors,
z-statistics, p- values, and the log-likelihood function value for each
equation.

##  5 Illustration

Markov chain models are used in interdisciplinary areas, such as economics,
business, biology, and engineering, with applications to predict long-term
behavior from traffic flow to stock market movements, among others. Modeling
and predicting stock markets returns is particularly relevant for investors
and policy makers. Since the stock market is a volatile environment, and the
returns are difficult to predict, estimating the set of probabilities that
describe these movements, might provide relevant input. Additionally,
incorporating the effect of key macroeconomic variables could provide a more
accurate picture of this specific environment. The following empirical
illustration aims to model stock returns of two indexes as a function of the
interest rate spread, specifically the 10-Year Treasury Constant Maturity
Minus 3-Month Treasury Constant Maturity.

The interest rate spread is a key macroeconomic variable and provides valuable
information regarding the economy state. Specifically, it has been used to
forecast recessions as in Estrella, Mishkin (1996), Dombrosky, Haubrich
(1996), Chauvet, Senyuz (2016), Tian, Shen (2019) and McMillan (2021).
Generically, short-term yields are lower than long-term yields when the
economy is in expansion. On the other hand, short-term yields are higher than
long-term yields when the economy is in recession. The difference between
these yields (or, more specifically, the yield curve’s slope) can be used to
forecast the state of the economy. Hence, this indicator might provide
relevant input for investors.

We considered the 5-week-day daily stock returns ($r_{t}=100\times
log(P_{t}/P_{t-1})$, where $P_{t}$ is the adjusted close price) of two
indexes, S&P500 and djia, from November $11^{th}$ 2011 to September $1^{st}$
2021 (2581 observations). Additionally, we considered the interest rate spread
($spread_{t}$), the 10-Year Treasury Constant Maturity Minus 3-Month Treasury
Constant Maturity. The data was retrieved from fred. Below in table 1, we have
the descriptive statistics of these variables.

Variable | Minimum | 1st Qu. | Median | Mean | 3rd Qu. | Maximum
---|---|---|---|---|---|---
$spread_{t}$ | -0.520 | 0.920 | 1.540 | 1.454 | 2.030 | 2.970
$r_{t;SP500}$ | -12.765 | -0.3198 | 0.070 | 0.054 | 0.518 | 8.968
$r_{t;DJIA}$ | -13.841 | -0.327 | 0.071 | 0.046 | 0.508 | 10.764
Table 1: Descriptive statistics

Moreover, to apply the model proposed, it is necessary to have a categorical
time series, thus we applied the following procedure:

$S_{st}=\begin{cases}1,r_{t}\leq\hat{q}_{s;0.25}\\\
2,\hat{q}_{s;0.25}<r_{t}<\hat{q}_{s;0.75}\\\ 3,r_{t}\geq\hat{q}_{s;0.75}\\\
\end{cases}$

where $\hat{q}_{s;\alpha}$ is the estimated quantile of order $\alpha$ of the
marginal distribution of $r_{t}$.

In Figures 5 and 6, we have the smoothed conditional probabilities of both
series. The number of observations is high, and the probabilities varied
abruptly in a small time frame, making the plots hard to read. So, to
simplify, a moving average model of order 5 (due to the frequency of the data)
was adjusted to these probabilities to illustrate how they evolve throughout
time. We see a similar behavior within each series regardless of whether it
depends on the previous states of $S_{1t}$ or $S_{2t}$. Additionally, the
scales of the graphs are small, indicating that these probabilities vary
around the same set of values.

(a) Probabilities of series 1 (SP500) depending on $spread_{t-1}$ and on
series 1 (SP500) previous state

(b) Probabilities of series 1 (SP500) depending on $spread_{t-1}$ and on
series 2 (DJIA) previous state

Figure 5: Conditional Probabilities of SP500’s series

(a) Probabilities of series 2 (DJIA) depending on $spread_{t-1}$ and on series
1 (SP500) previous state

(b) Probabilities of series 2 (DJIA) depending on $spread_{t-1}$ and on series
2 (DJIA) previous state

Figure 6: Conditional Probabilities of DJIA’s series

From this set of probabilities, we can estimate the model. Considering the
first equation, the effect of the probabilities depending on S&P500’s previous
state and the interest rate spread has a higher weight on the overall
probability. Also, this estimate is highly significant, presenting a $p$-value
close to zero. The effect of DJIA’s previous state in S&P500 is lower but it
is also significant for a 10% significance level. In the second equation, the
effect of S&P500’s previous state is higher than DJIA’s and both estimates are
highly significant.

    R> mmcx(cbind(s1, s2), bd$x_1, initial = c(1,1))

    $‘Equation 1‘
      Estimate Std. Error t value Pr(>|t|)
    1 0.685660   0.171241   4.004    0.000 ***
    2 0.314340   0.171241   1.836    0.066 *

    $‘LogLik 1‘
              [,1]
    [1,] -2636.355

    $‘Equation 2‘
      Estimate Std. Error t value Pr(>|t|)
    1 0.629992   0.176310   3.573    0.000 ***
    2 0.370008   0.176309   2.099    0.036 **

    $‘LogLik 2‘
              [,1]
    [1,] -2636.622

One of the advantages of this approach is the possibility to assess the
transition probabilities for specific values of $x_{t}$, in this case, the
interest rate spread.So, for both series, we calculated the transition
probabilities for this variable’s minimum and maximum value in the sample,
which are -0.52 and 2.97, respectively. In Figure 7, we have the transition
probabilities network for S&P500, corresponding to the minimum and maximum
value of the spread.

The most noticeable difference between these two networks is regarding the
transition probability from the second state to the third state. For the
maximum value of $spread_{t-1}$, the transition probability from the second
state to the third state is 0.6. So, when the economy is strong, one might
expect to have higher returns, when $t-1$ was in the second state. However,
this scenario shifts when considering the minimum value of $spread_{t-1}$. The
probability of obtaining higher returns, that is, being in state three,
becomes almost evenly distributed, regardless of the state in $t-1$. This
indicates the instability of the stock market, when the economy is weaker.
Another difference in these networks, is regarding the transition probability
from the third state to the first state. For the maximum value of
$spread_{t-1}$, this probability is 0.27 and for the minimum value increases
to 0.44. This is also expected, since when the economy is weaker, the
probability of having lower returns is greater.

(a) Conditional probabilities of series 1 for the maximum value of
$spread_{t-1}$

(b) Conditional probabilities of series 1 for the minimum value of
$spread_{t-1}$

Figure 7: Transition Probabilities of Series 1: S&P500

Considering the second equation, corresponding to the djia’s returns, we see a
similar behaviour as in S&P500’s networks. The transition probability from the
second state to the third state is higher for the maximum value of
$spread_{t-1}$ and the transition probability from the third state to the
first state is higher when we consider the minimum value of $spread_{t-1}$.
Although, the difference of this last probability between the minimum and
maximum value of $spread_{t-1}$ is not as big as in S&P500. Overall, the rest
of the probabilities structure, remains the same.

(a) Conditional probabilities of series 2 for the maximum value of
$spread_{t-1}$

(b) Conditional probabilities of series 2 for the minimum value of
$spread_{t-1}$

Figure 8: Transition Probabilities of Series 2: DJIA

##  6 Conclusions, limitations and further research

Several proposals for including of exogenous variables in MMC models have been
presented. The main limitations were associated with the high complexity of
the models to be developed and estimated. Additionally, most models considered
only categorical exogenous variables, existing a lack of focus on continuous
exogenous variables.

This work proposes a new approach to include continuous exogenous variables in
Ching et al. (2002) model for multivariate Markov chains. This is relevant
because it allows studying the effect of previous series and exogenous
variables on the transition probabilities.

The model is based on Ching et al. (2002) MMC model but considers non-
homogeneous Markov chains. Thus, the probabilities that compose the model are
dependent on exogenous variables. These probabilities are estimated as a usual
non-homogeneous Markov chain through a multinomial logit model. The model
parameters are then estimated through MLE, as well as the standard errors. We
developed a package with the estimation function of the model proposed. In
this, we considered the Augmented Lagrangian optimization method for
estimating the parameters through MLE. Additionally, we designed a Monte Carlo
simulation study to assess this model’s test power and dimension. The results
showed that the model detected a non-homogeneous Markov chain. Moreover, an
empirical illustration demonstrated the relevance of this new model by
estimating the probability transition matrix for different exogenous variable
values. Ignoring the effect of exogenous variables in MMC means that we would
not detect the probabilities’ changes according to the covariates’ values. In
this setting, one would have a limited view of the studied process. Hence,
this approach allows us to understand how a specific variable influences a
specific process.

The main contributions of this work are the development of a package with
functions for multivariate Markov chains, addressing the statistical inference
in these models and the inclusion of covariates. The limitations are related
to the implementation in R, specifically the optimization algorithm applied is
not common for MMC models, in that sense, it would be beneficial to study new
approaches to optimizing the maximum likelihood function as further research.
Additionally, extending this generalization to the MTD-probit model proposed
by Nicolau (2014) would also be relevant, which removes the constraints of the
model’s parameters and allows the model to detect negative effects.

## References

  * Adke, Deshmukh (1988) Adke S.R., Deshmukh S.R. Limit Distribution of a High Order Markov Chain  // Journal of the Royal Statistical Society. Series B (Methodological). 1988. 50, 1. 105–108.
  * Azzalini (1994) Azzalini A. Logistic regression for autocorrelated data with application to repeated measures  // Biometrika. 1994. 81, 4. 767–775.
  * Bartholomew (1968) Bartholomew J. Stochastic Models for Social Processes  // The Australian and New Zealand Journal of Sociology. 1968. 4, 2. 171–172.
  * Berchtold (1995) Berchtold A. Autoregressive Modelling of Markov Chains  // Proc. 10th International Workshop on Statistical Modelling. 1995. 104. 19–26.
  * Berchtold (1996) Berchtold A. Modélisation autorégressive des chaînes de Markov : utilisation d’une matrice différente pour chaque retard  // Revue de Statistique Appliquée. 1996. 44, 3. 5–25.
  * Berchtold (2001) Berchtold A. Estimation in the Mixture Transition Distribution Model  // Journal of Time Series Analysis. 2001. 22, 4. 379–397.
  * Berchtold (2003) Berchtold A. Mixture transition distribution (MTD) modeling of heteroscedastic time series  // Computational Statistics and Data Analysis. 2003. 41, 3-4. 399–411.
  * Berchtold et al. (2020) Berchtold A., Maitre O., Emery K. Optimization of the mixture transition distribution model using the march package for R  // Symmetry. 2020. 12, 12. 1–14.
  * Billingsley (1961) Billingsley P. The Lindeberg-Lévy Theorem for Martingales  // Proceedings of the American Mathematical Society. 1961. 12, 5. 788–792.
  * Bolano (2020) Bolano D. Handling covariates in markovian models with a mixture transition distribution based approach  // Symmetry. 2020. 12, 4.
  * Chauvet, Senyuz (2016) Chauvet M., Senyuz Z. A dynamic factor model of the yield curve components as a predictor of the economy  // International Journal of Forecasting. 2016\. 32, 2. 324–343.
  * Chen, Lio (2009) Chen D. G., Lio Y. L. A Novel Estimation Approach for Mixture Transition Distribution Model in High-Order Markov Chains  // Communications in Statistics - Simulation and Computation. 2009. 38, 5. 990–1003.
  * Ching et al. (2002) Ching W. K., Fung E. S., Ng M. K. A multivariate Markov chain model for categorical data sequences and its applications in demand predictions  // IMA Journal of Management Mathematics. 2002. 13, 3. 187–199.
  * Ching et al. (2003) Ching W. K., Fung E. S., Ng M. K. A Higher-Order Markov Model for the Newsboy’s Problem  // The Journal of the Operational Research Society. 2003. 54, 3. 291–298.
  * Ching, Ng (2006) Ching W. K., Ng M. K. Markov Chains: Models, Algorithms and Applications. 2006.
  * Ching et al. (2008) Ching W. K., Ng M. K., Fung E. S. Higher-order multivariate Markov chains and their applications  // Linear Algebra and its Applications. 2008. 428, 2-3. 492–507.
  * Damásio (2013) Damásio B. Multivariate Markov Chains - Estimation, Inference and Forecast. A New Approach: What If We Use Them As Stochastic Covariates? 2013.
  * Damásio, Mendonça (2019) Damásio B., Mendonça S. Modelling insurgent-incumbent dynamics: Vector autoregressions, multivariate Markov chains, and the nature of technological competition  // Applied Economics Letters. 2019. 26, 10. 843–849.
  * Damásio, Mendonça (2020) Damásio B., Mendonça S. Leader-follower dynamics in real historical time: A Markovian test of non-linear causality between sail and steam (co-)development, mimeo. 2020.
  * Damásio, Nicolau (2014) Damásio B., Nicolau J. Combining a regression model with a multivariate Markov chain in a forecasting problem  // Statistics & Probability Letters. 2014. 90. 108–113.
  * Damásio, Nicolau (2020) Damásio B., Nicolau J. Time inhomogeneous multivariate Markov chains : detecting and testing multiple structural breaks occurring at unknown dates. 2020.  (REM Working Papers).
  * Damásio et al. (2018) Damásio B., Louçã F., Nicolau J. The changing economic regimes and expected time to recover of the peripheral countries under the euro: A nonparametric approach  // Physica A: Statistical Mechanics and its Applications. 2018\. 507. 524–533.
  * Dombrosky, Haubrich (1996) Dombrosky A. M., Haubrich J. Predicting real growth using the yield curve  // Economic Review. 1996. I, Q. 26–35.
  * Estrella, Mishkin (1996) Estrella A., Mishkin F. S. The yield curve as a predictor of U.S. recessions  // Current Issues in Economics and Finance. 1996. 2, Jun.
  * Henningsen, Toomet (2011) Henningsen A., Toomet O. maxLik: A package for maximum likelihood estimation in R  // Computational Statistics. 2011. 26, 3. 443–458.
  * Islam et al. (2004) Islam M. A., Arabia S., Chowdhury R. I. A Three State Markov Model for Analyzing Covariate Dependence  // International Journal of Statistical Sciences. 2004. 3, i. 241–249.
  * Islam, Chowdhury (2006) Islam M. A., Chowdhury R. I. A higher order Markov model for analyzing covariate dependence  // Applied Mathematical Modelling. 2006. 30, 6. 477–488.
  * Jacobs, Lewis (1978) Jacobs P.A., Lewis A.W. Discrete Time Series Generated by Mixtures II : Asymptotic Properties  // Journal of the Royal Statistical Society: Series B (Methodological). 1978. 40, 2. 222–228.
  * Kalbfleisch, Lawless (1985) Kalbfleisch J. D., Lawless J. F. The analysis of panel data under a Markov assumption  // Journal of the American Statistical Association. 1985. 80, 392. 863–871.
  * Kijima et al. (2002) Kijima M., Komoribayashi K., Suzuki E. A multivariate Markov model for simulating correlated defaults  // Journal of Risk. 07 2002. 4.
  * Le et al. (1996) Le N. D., Martin R. D., Raftery A. Modeling Flat Stretches, Brusts, and Outliers in Time Series Using Mixture Transition Distribution Models  // Journal of the American Statistical Association. 1996. 91, 436. 1504–1515.
  * Logan (1981) Logan J. A structural model of the higher-order Markov process incorporating reversion effects  // The Journal of Mathematical Sociology. 1981. 8, 1. 75–89.
  * Lèbre, Bourguignon (2008) Lèbre S., Bourguignon P. Y. An EM algorithm for estimation in the mixture transition distribution model  // Journal of Statistical Computation and Simulation. 2008. 78, 8. 713–729.
  * Maitre, Emery (2020) Maitre O., Emery K. march: Markov Chains. 2020.  R package version 3.3.2.
  * Martin, Raftery (1987) Martin R. D., Raftery A. Non-Gaussian State-Space Modeling of Nonstationary Time Series: Comment: Robustness, Computation, and Non-Euclidean Models  // Journal of the American Statistical Association. 1987. 82, 400. 1044–1050.
  * McMillan (2021) McMillan D. G. Predicting GDP growth with stock and bond markets: Do they contain different information?  // International Journal of Finance & Economics. 2021. 26, 3. 3651–3675.
  * Mehran (1989) Mehran F. Analysis of Discrete Longitudinal Data: Infinite-Lag Markov Models  // Statistical Data Analysis and Inference. Amsterdam: North-Holland, 1989. 533–541.
  * Muenz, Rubinstein (1985) Muenz L. R., Rubinstein L. V. Markov Models for Covariate Dependence of Binary Sequences  // Biometrics. 1985. 41, 1. 91–101.
  * Newey, Mcfadden (1994) Newey W.K, Mcfadden D. Chapter 36 Large sample estimation and hypothesis testing  // Handbook of Econometrics. 4. 1994. 2111–2245.
  * Nicholson (2013) Nicholson W. DTMCPack: Suite of functions related to discrete-time discrete-state Markov Chains. 2013.  R package version 0.1-2.
  * Nicolau (2014) Nicolau J. A new model for multivariate markov chains  // Scandinavian Journal of Statistics. 2014. 41, 4. 1124–1135.
  * Nicolau, Riedlinger (2014) Nicolau J., Riedlinger F. I. Estimation and inference in multivariate Markov chains  // Statistical Papers. 2014. 56, 4. 1163–1173.
  * Pegram (1980) Pegram G. An Autoregressive Model for Multilag Markov Chains  // Journal of Applied Probability. 1980. 17, 2. 350–362.
  * Raftery (1985) Raftery A. A Model for High-Order Markov Chains  // Journal of the Royal Statistical Society: Series B (Methodological). 1985. 47, 3. 528–539.
  * Raftery, Tavaré (1994) Raftery A., Tavaré S. Estimation and Modelling Repeated Patterns in High Order Markov Chains with the Mixture Transition Distribution Model  // Applied Statistics. 1994. 43, 1. 179–199.
  * Rajarshi (2013) Rajarshi M.B. Statistical Inference for Discrete Time Stochastic Processes. 2013.
  * Regier (1968) Regier M. H. A Two-State Markov Model for Behavioral Change  // Journal of the American Statistical Association. 1968. 63, 323. 993–999.
  * Siu et al. (2005) Siu T. K., Ching W. K., Fung E. S., Ng M. K. On a multivariate Markov chain model for credit risk measurement  // Quantitative Finance. 2005. 5, 6. 543–556.
  * Spedicato (2017) Spedicato G. A. Discrete Time Markov Chains with R  // The R Journal. 07 2017.  R package version 0.6.9.7.
  * Spilerman, Singer (1976) Spilerman S., Singer B. The Representation of Social Processes by Markov Models  // American Journal of Sociology. 1976. 82, 1. 1–54.
  * Tian, Shen (2019) Tian R., Shen G. Predictive power of Markovian models: Evidence from US recession forecasting  // Journal of Forecasting. 2019. 38, 6. 525–551.
  * Venables, Ripley (2002) Venables W. N., Ripley B. D. Modern Applied Statistics with S. New York: Springer, 2002. Fourth.  ISBN 0-387-95457-0.
  * Wang et al. (2014) Wang C., Huang T. Z., Ching W. K. A new multivariate Markov chain model for adding a new categorical data sequence  // Mathematical Problems in Engineering. 2014\. 2014.
  * Wasserman (1980) Wasserman S. Analyzing social networks as stochastic processes  // Journal of the American Statistical Association. 1980. 75, 370. 280–294.
  * Wong, Li (2001) Wong C. S., Li W. K. On a mixture autoregressive conditional heteroscedastic model  // Journal of the American Statistical Association. 2001. 96, 455. 982–995.
  * Zhang et al. (2006) Zhang X., King M. L., Hyndman R. J. A Bayesian approach to bandwidth selection for multivariate kernel density estimation  // Computational Statistics and Data Analysis. 2006. 50, 11. 3009–3031.
  * Zhu, Ching (2010) Zhu D. M., Ching W. K. A new estimation method for multivariate Markov chain model with application in demand predictions  // Proceedings - 3rd International Conference on Business Intelligence and Financial Engineering, BIFE 2010. 2010. 126–130.


# Can Noise on Qubits Be Learned in Quantum Neural Network? A Case Study on
QuantumFlow

Zhiding Liang†, Zhepeng Wang‡, Junhuan Yang¶, Lei Yang¶, Jinjun Xiong§, Yiyu
Shi†, Weiwen Jiang‡ †University of Notre Dame, IN, USA. ‡George Mason
University, VA, USA.
¶University of New Mexico, NM, USA. §University at Buffalo, NY, USA.
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In the noisy intermediate-scale quantum (NISQ) era, one of the key questions
is how to deal with the high noise level existing in physical quantum bits
(qubits). Quantum error correction is promising but requires an extensive
number (e.g., over 1,000) of physical qubits to create one “perfect” qubit,
exceeding the capacity of the existing quantum computers. This paper aims to
tackle the noise issue from another angle: instead of creating perfect qubits
for general quantum algorithms, we investigate the potential to mitigate the
noise issue for dedicate algorithms. Specifically, this paper targets quantum
neural network (QNN), and proposes to learn the errors in the training phase,
so that the identified QNN model can be resilient to noise. As a result, the
implementation of QNN needs no or a small number of additional physical
qubits, which is more realistic for the near-term quantum computers. To
achieve this goal, an application-specific compiler is essential: on the one
hand, the error cannot be learned if the mapping from logical qubits to
physical qubits exists randomness; on the other hand, the compiler needs to be
efficient so that the lengthy training procedure can be completed in a
reasonable time. In this paper, we utilize the recent QNN framework,
QuantumFlow, as a case study. Experimental results show that the proposed
approach can optimize QNN models for different errors in qubits, achieving up
to 28% accuracy improvement compared with the model obtained by the error-
agnostic training.

###### Index Terms:

Quantum machine learning, qubit mapping, noisy intermediate-scale quantum,
model training

##  I Introduction

Quantum computing has the quantum potential to provide exponential speedup
over the classical computers [8, 27, 1, 22, 2], which is considered to be one
of the best candidates for solving some of the complicated computational
problems that cannot be solved by classical computing in a reasonable time [8,
30, 9]. We are now witnessing the rapid development in both quantum hardware
and quantum algorithms: on the hardware side, the newly emerging quantum
computers (such as IBM Q [11]) are developed with superconductor and can
perform superposition and quantum entanglement on the physical quantum bits
(qubits); on the other side, different quantum algorithms have been devised,
such as Shor’s algorithm for primary factorization [23] and the quantum
approximate optimization algorithm (QAOA) for optimization problems [6].
However, the current quantum computers have a high noise level and
intermediate scale, known as noisy intermediate-scale quantum (NISQ). How to
deal with the noise with limited number of qubits becomes an imminent problem.

Quantum Error Correction (QEC) [4] is a promising solution to solve the noise
issue in quantum computing by creating “perfect qubits” (e.g., error rate to
be less than $10^{-15}$). It applies a set of noisy physical qubits to
represent one perfect qubit, which can correct the errors made by the noise
and significantly reduce the error rate. Once the QEC is achieved, it can be
beneficial for general quantum algorithms. But, there is still a long way to
achieve QEC, because it would require a large number (e.g., larger than 1,000
[20]) of physical qubits for a perfect qubit. On the other hand, quantum
computer in the NISQ era has a very limited number of qubits; for example, IBM
Q has at most 65 qubits on the IBM-Manhattan backend. Therefore, it is
critical to find a more efficient way to solve the noise issue in the NISQ
era.

Instead of creating perfect qubits for general quantum algorithms with QEC, in
this paper, we aim to address the noise issue in an application-specific way,
such that no more qubits or a very small number of additional qubits are
needed, while the effects of error can be mitigated. More specifically, we
would like to explore the design of a dedicated quantum algorithm to make it
resilient to noise. Quantum neural network (QNN) is a good candidate to
achieve such a goal since the neural network in computing-in-memory (CiM)
platforms has already demonstrated its capability to be resilient to noise
[16]. Compared with the design of CiM-based neural network, quantum computing
brings more challenges to identify the robust QNNs: (1) in terms of different
properties (e.g., topology, error rate, number of qubits) on quantum
computers, the implementations (or compilation) of the same QNN on hardware
can be different, making the learned model useless; (2) the mapping of QNN to
quantum computer needs to be integrated into the training process, leading to
lengthy training time; and (3) the update of one weight in the noisy quantum
computer may affect other weights, resulting in the traditional training
approach to be not applicable. The network structure can affects robust [17,
31, 32, 15, 14, 24, 29, 3, 28], which we will explore in the future.

In this paper, we base on a recent QNN design, i.e., QuantumFlow [12, 13], to
demonstrate how to overcome the above challenges. We first devise an
application-specific mapping algorithm. The proposed approach can fix the
mapping from logical circuits to physical qubits for different quantum
computers; meanwhile, it can minimize the number of gates to be used in
physical qubits. In addition, the application-specific mapping follows the
same procedure, which can significantly reduce the optimization time to better
support training the neural network. Then, based on the proposed mapping
algorithm, we build up an error-aware training framework to identify the
robust QNN model. Kindly note that the training procedure is completed at
classical computing in our implementation, which can be integrated into
quantum computing when the scale of the quantum computer is ready for QNNs.

The contribution of this paper is three-fold as follows.

  * •

To our best knowledge, this is the very first work to learn the quantum error
in quantum neural networks, and we demonstrate that a robust quantum neural
network, namely QF-RobustNN, can be obtained.

  * •

To support learning qubits’ noise in quantum neural networks, we first devise
an application-specific mapping method for one state-of-the-art quantum neural
network, QuantumFlow.

  * •

Then, a general training framework is proposed to conduct the error-aware
learning, which can be easily extended to support training on the quantum
computer in the future when the quantum computer is ready for quantum neural
networks, in terms of the number of qubits.

The paper is organized as follows: In Section II, we present the preliminaries
and motivations of this work. Section III presents our proposed mapping
algorithm and training framework. Evaluation results are reported in Section
IV. Section V concludes this paper.

Figure 1: BRIDGE and SWAP: (a) A CNOT gate with distance 2. (b) Assumed
physical qubits connection. (c) Use the BRIDGE gate to achieve gate (a). (d)
Use SWAP gate to achieve gate (a).

##  II Preliminaries

Qubits: In classical computing, the basic unit to represent data is a binary
bit (i.e., 0 or 1); similarly quantum computing utilizes quantum bit (qubit)
for data representation. Unlike a classical bit that can only represent 0 or
1, the qubit is to represent 0 and 1 simultaneously. More specifically, it
describes a state $|\psi\rangle$ with the combination of the 0-state (i.e.,
$|0\rangle$) and 1-state (i.e., $|1\rangle$); i.e.,
$|\psi\rangle=\alpha|0\rangle+\beta|1\rangle$, where $\alpha$ and $\beta$ are
the amplitudes of 0-state and 1-state, and $\alpha^{2}$ and $\beta^{2}$
represent their probabilities. When $\alpha\neq 0$ and $\beta\neq 0$, the
qubit can represent two-state at the same time, known as the “superposition”
state.

Quantum Gates: In analogy to the classical logic gates that process classical
data, quantum gates manipulate quantum information, i.e., the quantum states
qubits. A quantum system is typically initialized to zero state (i.e.,
$|\psi\rangle=|0\rangle$). Given a Hamiltonian, [7] has demonstrated that the
Hamiltonian can be solved by Schrödinger’s equation, which can be decomposed
into various combinations of single-qubit (e.g., H, S, and T gates) and multi-
qubit quantum gates (e.g., CNOT gate). It has been proven that the combination
of Clifford gates (e.g., CNOT, H, S gates) and T gate is a universal gate set
to represent arbitrarily complex gates for quantum computing.

Quantum Circuit: A series of logical gates can construct the classical
circuit; likewise, a series of quantum gates form the quantum circuit.
According to the function/algorithm needing to be implemented, a quantum
circuit can be is formed by using the quantum gates in the universal gate set.
The length/depth of the circuit reflects the time complexity of the quantum
version of the implemented algorithm, and the width of the circuit indicates
the number of qubits to be used.

Quantum Processor: The quantum processor (a.k.a., quantum computer) is built
upon a set of physical qubits, which are implemented by superconductor [18],
IoN-Trap [19], etc. Taking the commonly superconductor-based quantum processor
as an example, it has two limitations: (1) qubits are not fully connected; and
(2) noise exists in operating qubits or a pair of qubits.

TABLE I: Accuracy in Different Noise Models Noise Model | Affected Gates | Error Rate | Accuracy | Time
---|---|---|---|---
No Error | 0 | - | 0.9804 | 5.00s
Flip Error | X, CX, CCX | 0.1 | 0.5333 | 568.50s
Flip Error | X, CX, CCX | 0.01 | 0.8824 | 540.00s
Phase Error | Z, CZ | 0.1 | 0.6429 | 545.00s
Phase Error | Z, CZ | 0.01 | 0.9167 | 511.14s
Flip+Phase | X,CX,CCX,Z,CZ | 0.1 | 0.451 | 628.00s
Flip+Phase | X,CX,CCX,Z,CZ | 0.01 | 0.7778 | 532.04s

SWAP and BRIDGE for Limitation (1). Figure 1 shows an example with 3 qubits
organized in a linear form, where logical qubit $|q_{i}\rangle$ is mapped to
physical qubit $p_{i}$, i.e., $M(|q_{i}\rangle)=p_{i}$. And a CNOT gate
operates between $M(|q_{0}\rangle)$ and $M(|q_{2}\rangle)$. Obviously, we
cannot directly conduct CNOT the physical qubits has no channel for
communication. To solve this, a BRIDGE operation and a SWAP operation in
Figure 1(c)-(d) are designed. The idea of BRIDGE in Figure 1(c) is to create a
unitary operator that performs the function of CNOT gate between two
physically separate qubits, while the state of the middle qubit (i.e.,
$M(|q_{1}\rangle)$) is not changed. On the other hand, the SWAP in Figure 1(d)
is to exchange the quantum states between $p_{0}$ and $p_{1}$, as such, after
the SWAP, we will have $M(|q_{0}\rangle)=p_{1}$ and $M(|q_{1}\rangle)=p_{0}$,
which can perform CNOT on $|q_{0}\rangle$ and $|q_{1}\rangle$. It is worth
noticing that both BRIDGE and SWAP can work for the communication between two
qubits needing to cross more than 1 physical qubit, but needing more CNOT
gates for implementation, introducing higher cost.

Noisy Qubits leading to Limitation (2): Due to the physical material, there
exists a high noise level in quantum qubits. In classical computing, the “bit
flip” is a typical error; when it comes to the qubit, there are more types of
error, which can be concluded as follows:

(1) Qubit-flip Error: This is manifested as a state flip of a qubit from
$|0\rangle$ to $|1\rangle$ or versus.

(2) Phase Error: This is manifested as a phase shift of a qubit, which can be
regarded as the flip on the phase (e.g., perform a Z gate).

(3) Readout Error: This is specified by a list assignment probability vector
$P(n|m)$, assuming $m$ is a bit string from an ideal measurement, while $n$ is
from a noisy measurement. $P(n|m)$ reflects the number of errors exists in
$n$.

(4) Depolarizing Error: An n-qubit channel stores depolarization error
parameterized by a depolarization probability $p$.

Quantum Compiler: In order to implement a quantum circuit to a specific
quantum process, a compiler is essential to map the logical circuit to
physical qubits with the different objectives (e.g., minimizing the number of
quantum gates, circuit depth, fidelity, etc.). The existing quantum compiler
is almost designed for general purpose [5], but we observe that they may be
too costly to be applied for the scenario needing to compile a large number of
logical circuits. And this work aims to provide an application-specific
quantum compiler to find the most efficient mapping.

Target Quantum Application: The application considered in this work is Quantum
Neural Network (QNN), which is recently emerged and attacked wide attention
from both academics and industry. In the training phase of QNN, it may require
compiling an extensive number of quantum circuits, which has a high demand on
the compiler efficiency. For example, FFNN [25] and QuantumFlow [12] will
build different quantum circuits for different weights.

Figure 2: Illustration of error unpredictable caused by compilation: (a) and
(c) shows two logical quantum circuits L1 and L2. (b) and (d) describe the
steps processing in the physical circuit P1 and P2, which corresponding to L1
and L2, respectively.

##  III Motivation

Noise can significantly affect the performance of QNN. In Table I, we report
the results on the MINST database: QuantumFlow on perfect qubits can obtain a
98.04% of accuracy. When we add bit-flip errors with error rates of 0.1 and
0.01, the accuracy degraded to 53.33% and 88.24%, respectively; when we add
phase errors with error rates of 0.1 and 0.01, the accuracy is 64.29% and
91.67%, respectively Moreover, for a quantum processor with both bit-flip and
phase-error and the error rate to be 0.1 and 0.01, the accuracy further
reduced to 45.10% and 77.78%, respectively. From the above results, it is easy
to see how noise interferes with QuantumFlow. With such observation, we are
thinking of whether it is possible to learn the error in the training phase,
such that we can obtain a more robust QNN.

We have seen the success of learning errors in classical computing [10]. But,
when it comes to quantum computing, more challenges pose up.

Challenge 1: Error unpredictable on the quantum circuit.

In quantum processors, the error on qubits is usually different; if a weight
in QNN cannot be mapped to a fixed qubit for execution at run-time, it will
lead to the error being unpredictable. Consider a QNN with two different
weights, which are embedded into the circuit with different quantum circuits
(L1 in Figure 2(a) and L2 in 2(c)).

For L1, Figure 2(b) shows its implementation on the quantum process in Figure
1(b). Assume the logical qubits ($q_{0},q_{1},q_{2},q_{3}$) are mapped to the
physical qubits as ($p_{0},p_{1},p_{3},p_{4}$) after the data pre-processing
(blue box). In embedding the weights into the circuit (last CZ gate with two
dots), qubit $q_{0}$ is on physical qubit $p_{2}$, while $q_{1}$ is on
$p_{1}$. On the other hand, due to the change of weights, L2 does not have the
second CZ gate (blue dashed box). As a result, qubit $q_{0}$ is on physical
qubit $p_{1}$, while $q_{1}$ is on $p_{0}$. As a result, when we execute the
circuit, the error for the same weight can be totally different in two QNN
models. This will hinder us to learn the error in the QNN model.

TABLE II: Compilation time of QUEKO [26]. Quantum Gate(with 4 qubits) | Compilation Time
---|---
CCX | 1.73s
CZ | 0.14s
CCX + CZ +CCX | 16.43s
CCX + CZ +CCX + CCX | 91.91s
CCX + CZ +CCX + CCX + CZ | 128.20s
CCX + CZ +CCX + CCX + CZ + CCX | 439.67s

Challenge 2: Existing compiler may not for training.

It is well known that compiling is an NP-Hard problem. For a general-purpose
compiler, it will have a high execution time for better exploring the search
space. On the other hand, since we need to integrate the compiling into the
training iteration, and it is also well known that training in the neural
network is extremely time-consuming. In terms of compilation time, we exam one
of the state-of-the-art compilers [26] with the configuration of “fidelity”
and “transition” (a high-speed processing mode). Table II reports the result:
the compilation time increases significantly with the increase of quantum gate
numbers in the quantum circuit. More specifically, it takes less than 1 second
for one gate, but the compilation time grows to more than 400 seconds for 6
gates. This is commonly existing different general-purpose compiler since the
increase of the number of gates can exponentially increase the search space,
leading the compilation time increase exponentially. However, the quantum
circuit in QuantumFlow is far more complicated than the circuit with 6 gates
in Table II, and we need to conduct such compilation each time for a new
weight is identified in the training phase.

Motivation: After summarizing the above-described observations and challenges,
it is clear that we need to resolve the effects of qubits’ error to make QNN
effective, and we need to reduce the compilation time to make it possible to
learn the error in the training phase.

##  IV Proposed Algorithm

###  IV-A Problem Statement

With the demonstration of the existing challenges and our observations, it is
clear that the existing errors in quantum computers will greatly reduce
accuracy. In addition, in order to learn the error in the QNN mode, the
compilation should be more efficient to support the training. We formally
define the problem of ‘learning quantum error in QNN’ as follow:

Figure 3: Illustration of the proposed training framework to generate QF-
RobustNN.

Quantum-error-aware neural network training: Given a neural network and its
implementation on physical qubits, how to search for the weights of the neural
network, such that the accuracy of the neural network model can be maximized
with the consideration of error on qubits.

In the above problem, the implementation of a given quantum circuit needing to
be implemented by a compiler, and therefore, we will further study how to
design an application-specific compiler to support the training procedure,
which is defined as follows.

Application-specific compiler: Given a dataset, the quantum algorithm, and
quantum processor, how to map logical circuit to the quantum processor, such
that each qubit is fixed to be mapped to a physical qubit and the circuit
length can be minimized.

###  IV-B Quantum-error-aware training framework

In this section, we provide detailed procedures for the operation of the
quantum-error-aware training framework, as shown in Figure 3. The framework is
composed of four components: i) the start point with trained weights; ii) the
logical quantum circuit created by QuantumFlow; iii) the mapping from logical
circuit to physic qubits; and iv) complete the inference to obtain model
accuracy and tune the parameters accordingly. If the model is converged after
the training step, we will terminate the algorithm and output the identified
QF-RobustNN.

To better demonstrate the framework, we formulate the procedures as follows:

$\displaystyle c_{i}=circ(W_{i})$ (1) $\displaystyle m_{i}=Map(c_{i},PhyQ)$
$\displaystyle e_{i}=Error(m_{i})$ $\displaystyle
a_{i}=Inference(m_{i},e_{i})$

For one identified weight in the $i^{th}$iteration, the circuit $c_{i}$ is
obtained, where the $circ$ represents obtain the circuit from a QNN framework
(e.g., QuantumFlow). Then, we map the circuit to the physical qubits (PhyQ)
with a mapping function $Map$, which will be introduced in the next
subsection. After this, we will build up the error model $e_{i}$ based on the
mapping $m_{i}$. Finally, we will conduct the inference to obtain the accuracy
($a_{i}$) of the QNN with weight $W_{i}$, which is obtained by executing the
physical circuit $m_{i}$ with the error model $e_{i}$. After this, we will be
going through the training procedure to update the weights in the QNN.

➀ Train QNN to learn error info: In the third step in Figure 3, the QNN
training model is built up with the purpose of searching for the best training
weights is created. The training weight will be learned with the error
function $Error$ to become be aware of the quantum error and make the model
resilient to the error $Error$. Given an initial weight that performs best in
perfect quantum qubits, the training procedure is to update the weights. Here,
we can apply different optimization approach to update the weights. For
example, we can apply deep reinforcement learning, where a recurrent neural
network (RNN)-based controller will guide the updates. More specifically, the
accuracy $a_{i}$ will be the input of the controller to update the parameters
of the RNN, and the RNN will then generate the weights to update the training
weights of the first block in Figure 3. The best performing weight is denoted
as ’Searched_Weight’ by choosing the highest accuracy among all results.
Comparing the results of ’Searched_Weight’ with the results of the initial
weights (called baseline weights), we can evaluate the effectiveness of the
proposed training framework, which will be reported in Section V.

Figure 4: Describe the proposed qubits mapping algorithm steps and the
corresponding gate connective behavior. (a) CCCZ gate and the decomposed
version of CCCZ gate. (b) CCX gate and the decomposed version of CCX gate. (c)
is the assumed physical qubit coupling graph and the step-by-step operations
of the proposed qubit mapping algorithm.

➁ Application-specific Mapping $C^{N}Z$ gates are the common gates utilized by
QuantumFlow to implement the quantum neural networks, where $N$ represents
there are $N$ control end with one target qubits using $Z$ gate. To simplify
the presentation, we choose $C^{3}Z$ (i.e., CCCZ) gates as an example to
describe our application-specific qubit mapping algorithm.

The CCCZ gate and the decomposed circuit of the CCCZ gate are illustrated in
Figure 4(a). Then, the decomposed circuit of the CCX gate which shows in
Figure 4(b) that consists of basic CNOT gates to be executed during CCX gate
operation that revealing logic qubits connection relationship. The topology of
the physical qubits and logical-to-physical qubits initial mapping is provided
in Figure 4(c). The assumed distribution and connectivity of physical qubits
are shown in the coupling graph, where $Q_{n}$ refers to the logic qubits
$q_{n}$ corresponding to their physical qubits in the quantum circuit and
$AUX_{n}$ refers to the auxiliary qubits $aux_{n}$ corresponding to their
physical qubits in a quantum circuit.

To make the qubit mapping for different $C^{N}Z$ gates to be fixed, our design
philosophy is to make the logical-to-physical mapping to be the same at the
beginning and end of the mapping procedure, meanwhile, we do not incur
additional SWAP gates to minimize the circuit length. The procedure of the
mapping algorithm is shown in Figure 4(d). The dotted boxes with the same
color in Figure 4(a) and Figure 4(d) have correspondence. In the first step,
we need first to conduct the qubit level logical-to-physical mapping. This
will be followed by the rule of interleaving the computing qubits and
auxiliary qubits, and they are mapped to a chain of qubits. For example,
Figure 4(c) gives an initial logical-to-physical mapping:
$\\{q_{0},aux_{0},q_{1},aux_{1},q_{2},q_{3}\\}\rightarrow\\{P_{0},P_{1},P_{2},P_{3},P_{4},P_{5}\\}$.
For convenient corresponding understanding, we refer
$\\{P_{0},P_{1},P_{2},P_{3},P_{4},P_{5}\\}$ as
$\\{Q_{0},AUX_{0},Q_{1},AUX_{1},Q_{2},Q_{3}\\}$ in Figure 4(d). Then, the
computation of CCX gates will be arranged in two-phase: forward phase and
backward phase.

In forward phase, at $t_{1}$, the CCX gate in black dotted box act on $q_{0}$,
$q_{1}$ and $aux_{0}$ in Figure 4(a). The corresponding decomposed circuit in
Figure 4(b) should be executed. It is obvious to see in the coupling graph in
Figure 4(c) that $Q_{0}$ has already been connected to $AUX_{0}$, and
$AUX_{0}$ is already connected to $Q_{1}$ by an edge of the coupling graph in
black dotted box, so that $g_{1}$, $g_{2}$, $g_{3}$, and $g_{4}$ in Figure
4(b) can be directly executed. However, we still need to build connection
between $Q_{0}$ and $Q_{1}$. Thus, we insert a SWAP operation
SWAP$(Q_{1},AUX_{0})$ before $g_{5}$ to move logic qubit $q_{1}$ toward
$aux_{0}$. After the SWAP operation, the mapping is update to
$\\{q_{0},aux_{0},q_{1},aux_{1},q_{2},q_{3}\\}\rightarrow\\{Q_{0},Q_{1},AUX_{0},AUX_{1},Q_{2},Q_{3}\\}$,
and $Q_{0}$ is connected to $Q_{1}$ that $g_{5}$ and $g_{6}$ can be executed.
At $t_{2}$, a same procedure as $t_{1}$ is performed, and finally we can see
that $|Q_{1}\rangle$ and $|AUX_{1}\rangle$ have connection in physical qubits,
and we can perform the CZ gate in place at $t_{3}$.

In backward phase, operations are processed at $t_{4}$ and $t_{5}$ in Figure
4(c). At $t_{4}$, The required connection relationship is $Q_{2}$ connected to
$AUX_{1}$, $AUX_{0}$ connected to $AUX_{1}$, and $Q_{2}$ connected to
$AUX_{0}$. First step is put a SWAP operation SWAP$(Q_{2},AUX_{1})$ in yellow
dotted box before executing $g_{1}$ to move logic qubit $q_{2}$ toward
$aux_{1}$. And the mapping is update to
$\\{q_{0},aux_{0},q_{1},aux_{1},q_{2},q_{3}\\}\rightarrow\\{Q_{0},Q_{1},AUX_{0},AUX_{1},Q_{2},Q_{3}\\}$.
This step is in order to build connection between $AUX_{0}$ and $AUX_{1}$,
$Q_{2}$ and $AUX_{1}$ by the edge of the coupling graph in yellow dotted box
in Figure 4(c), so that $g_{1}$, $g_{2}$, $g_{3}$ and $g_{4}$ can be executed.
The $AUX_{0}$ is needed to connect with $Q_{2}$. Here we do a BRIDGE operation
BRIDGE$(Q_{2},AUX_{0})$ in order to save the amounts of basic CNOT gates,
because if we do twice SWAP operation here, the cost number of basic CNOT
gates will be 8. BRIDGE operation could be only using 3 basic gate here due to
the middle logic qubit $|aux_{1}\rangle$ could be in state $|0\rangle$. And
after BRIDGE operation, $AUX_{0}$ is connected with $Q_{2}$ via the edge of
bridge, $g_{5}$ and $g_{6}$ in yellow dotted box can be executed. At $t_{5}$,
the same procedure is repeated. As a result, the qubit level logical-to-
physical mapping is resumed, so that we can achieve our goal to make all
weights to be mapped in a fixed and predictable way.

##  V Results and Discussions

This section describes the experimental setup and results. We are given an
initial weight that performs best on a synthetic dataset in the perfect model
as our baseline for accuracy comparison. In addition, the existing compilers
will be the baselines for efficiency comparison.

TABLE III: Accuracy Test Result From Proposed Algorithm in Different Models Error Rate | Baseline | QF-RobustNN
---|---|---
Acc. | Acc. | Weight
0 (Perfect) | 94% | - | [-1,-1,1,1,1,1,1,1], [1,1,1,1,1,1,1,1]
0.0001 | 84% | 84% | [-1,-1,1,1,1,1,1,1], [1,1,1,1,1,1,1,1]
0.0005 | 75% | 76% | [1,1,1,-1,-1,-1,-1,-1], [-1,-1,-1,1,1,-1,-1,-1]
0.001 | 74% | 76% | [1,1,-1,1,1,-1,-1,1], [-1,1,1,1,-1,1,1,1]
0.01 | 75% | 80% | [-1,-1,1,-1,1,-1,1,1], [1,1,1,-1,-1,1,-1,1]
0.05 | 49% | 75% | [1,1,-1,-1,1,-1,1,1], [-1,1,-1,-1,-1,-1,-1,1]
0.1 | 47% | 75% | [1,1,-1,-1,1,-1,1,-1], [1,1,-1,-1,1,-1,1,-1]
TABLE IV: Compilers Elapsed Time Result Compiler Name | Circuit Complexity Level
---|---
Simple | Middle | Complex
QUEKO [26] | 484.755s | 5332.903s | Over 10 hours
HA [21] | 4.765s | 4.696s | 4.201s
OURS | 0.008s | 0.015s | 0.020s
Imp. (vs QUEKO) | 60594.38$\times$ | 355526.87$\times$ | >1800000$\times$
Imp. (vs HA) | 595.63$\times$ | 313.07$\times$ | 210.05$\times$

A. Accuracy Comparison

We evaluate the performance of QF-RobustNN on both IBM Qiskit based simulator
and the IBM Quantum Processor (i.e., the ibmq_montreal backend with 27
qubits).

For the simulation, we set up different error rates with both bit-flip and
phase-shift for the evaluation. Table III reports the results. When the error
rate is 0.0001, the obtained QF-RobustNN is exactly the same as the baseline
model, because the baseline model is the one with the highest accuracy for
perfect qubits, and when the error rate is small, it still performs better
than other models. However, when the error rate increased to larger than
0.0005, we observe that QF-RobustNN becomes different from the baseline model.
Specifically, QF-RobustNN achieves 1%, 2%, and 5% accuracy improvements for
the error rates set to 0.0005, 0.001, 0.01, respectively. Furthermore, when
the error increased to 0.05 and 0.1, the baseline model has less than 50%
accuracy, but QF-RobustNN identified by the proposed error-aware training
framework can achieve 75% accuracy, with 24% and 28% accuracy gain. From the
above results, it is clear that the proposed training framework can work for
different error settings. In addition, with the increase of qubits’ error, the
model that works best for the perfect setting can easily become useless. This
emphasizes the importance of conducting error-aware learning.

B. Efficiency Comparison

We apply two existing compilers to be the baseline for efficiency comparison,
including (1) Quantum Mapping Examples with Known Optimal (QUEKO) [26] and (2)
heuristic-based Hardware-Aware mapping algorithm (HA) [21]. Table IV reports
the comparison results on elapsed time. Specifically, we run 3 different
circuits with circuit complex levels: simple, middle, and complex. QUEKO aims
to explore the best solution for the general-purpose quantum circuit, and it
takes more time than the other two compilers. Specifically, QUEKO takes
484.775 seconds and 5332.903 seconds for simple and middle circuits, while it
cannot optimize the complex circuit within 10 hours. On the other hand, HA can
complete all three circuits within 5 seconds. Furthermore, the proposed
application-specific compiler can complete the mapping within 0.1 seconds,
achieving over 200$\times$ speedup over HA. Results demonstrate that the
proposed application-specific compiler can be efficiently integrated into the
training framework.

C. Accuracy, Efficiency, and Circuit Depth

Finally, we execute benchmarks on ibmq_montreal to have an end-to-end
comparison between the solution obtained by using HA and QF-RobustNN obtained
by the proposed approaches. Table V reports the results. By using the baseline
model and our compiling approach, the resultant system provides the accuracy
of 74%. By applying QF-RobustNN, we can achieve improvement on accuracy to 75%
using HA and 80% using our application-specific compiler. One reason that we
can obtain higher accuracy is because our approach can better optimize the
circuit with less number of extra SWAP gates, where HA uses 43 SWAP gates
while the figure is 12 in our approach. In addition, we can see that our
proposed compiler can achieve better efficiency over HA.

TABLE V: Compilers Results on ibmq_montreal Compiler Name | Model | Accuracy | Extra SWAP gate | Elapsed Time
---|---|---|---|---
OURS | baseline | 74% | 15 | 0.020s
HA | QF-RobustNN | 75% | 43 | 3.955s
OURS | QF-RobustNN | 80% | 12 | 0.014s

##  VI Conclusion

In the NISQ era, noise is one critical issue needing to be resolved in quantum
computing in order to make quantum computing practical. It is still a long way
to go before we can achieve perfect qubits, either achieving a breakthrough in
material or realize quantum error correction. In such a background, we
proposed to design an application-specific optimization procedure to mitigate
the effects from the qubits’ noise. Targeting Quantum Neural Network, we
proposed the optimization framework to learn the error in the training
procedure. We further proposed an application-specific compiler to match the
needs of high efficient mapping to support the lengthy training procedure.
Results on both simulation and actual quantum computers demonstrate that such
an error-aware optimization can improve performance.

## Acknowledgment

Access to the IBM Q Network was obtained through the IBM Q Hub at NC State.
Special thanks to Zhirui Hu for her extensive discussion and help with this
work.

## References

  * [1] Frank Arute et al.  “Quantum supremacy using a programmable superconducting processor”  In _Nature_ 574.7779 Nature Publishing Group, 2019, pp. 505–510
  * [2] Koen Bertels et al.  “Quantum Accelerator Stack: A Research Roadmap”  In _arXiv preprint arXiv:2102.02035_ , 2021
  * [3] Jinglei Cheng, Haoqing Deng and Xuehai Qia  “Accqoc: Accelerating quantum optimal control based pulse generation”  In _2020 ACM/IEEE 47th Annual International Symposium on Computer Architecture (ISCA)_ , 2020, pp. 543–555  IEEE
  * [4] John Chiaverini et al.  “Realization of quantum error correction”  In _Nature_ 432.7017 Nature Publishing Group, 2004, pp. 602–605
  * [5] Alexander Cowtan et al.  “On the qubit routing problem”  In _arXiv preprint arXiv:1902.08091_ , 2019
  * [6] Edward Farhi, Jeffrey Goldstone and Sam Gutmann  “A quantum approximate optimization algorithm”  In _arXiv preprint arXiv:1411.4028_ , 2014
  * [7] David J Griffiths and Darrell F Schroeter  “Introduction to quantum mechanics”  Cambridge University Press, 2018
  * [8] Lov K Grover  “A fast quantum mechanical algorithm for database search”  In _Proceedings of the twenty-eighth annual ACM symposium on Theory of computing_ , 1996, pp. 212–219
  * [9] Aram W Harrow, Avinatan Hassidim and Seth Lloyd  “Quantum algorithm for linear systems of equations”  In _Physical review letters_ 103.15 APS, 2009, pp. 150502
  * [10] Vojtěch Havlı́ček et al.  “Supervised learning with quantum-enhanced feature spaces”  In _Nature_ 567.7747 Nature Publishing Group, 2019, pp. 209–212
  * [11] IBM  “IBM Q Experience Device”  In _https://www.reasearch.ibm.com/ ibm-q/technology/experience/_ , 2019
  * [12] Weiwen Jiang, Jinjun Xiong and Yiyu Shi  “A co-design framework of neural networks and quantum circuits towards quantum advantage”  In _Nature communications_ 12.1 Nature Publishing Group, 2021, pp. 1–13
  * [13] Weiwen Jiang, Jinjun Xiong and Yiyu Shi  “When Machine Learning Meets Quantum Computers: A Case Study”  In _2021 26th Asia and South Pacific Design Automation Conference (ASP-DAC)_ , 2021, pp. 593–598  IEEE
  * [14] Weiwen Jiang et al.  “Accuracy vs. efficiency: Achieving both through fpga-implementation aware neural architecture search”  In _Proceedings of the 56th Annual Design Automation Conference 2019_ , 2019, pp. 1–6
  * [15] Weiwen Jiang et al.  “Achieving super-linear speedup across multi-fpga for real-time dnn inference”  In _ACM Transactions on Embedded Computing Systems (TECS)_ 18.5s ACM New York, NY, USA, 2019, pp. 1–23
  * [16] Weiwen Jiang et al.  “Device-circuit-architecture co-exploration for computing-in-memory neural accelerators”  In _IEEE Transactions on Computers_ 70.4 IEEE, 2020, pp. 595–605
  * [17] Weiwen Jiang et al.  “Standing on the shoulders of giants: Hardware and neural architecture co-search with hot start”  In _IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems_ 39.11 IEEE, 2020, pp. 4154–4165
  * [18] C Monroe et al.  “Large-scale modular quantum-computer architecture with atomic memory and photonic interconnects”  In _Physical Review A_ 89.2 APS, 2014, pp. 022317
  * [19] Christopher Monroe and Jungsang Kim  “Scaling the ion trap quantum processor”  In _Science_ 339.6124 American Association for the Advancement of Science, 2013, pp. 1164–1169
  * [20] Prakash Murali et al.  “Noise-adaptive compiler mappings for noisy intermediate-scale quantum computers”  In _Proceedings of the Twenty-Fourth International Conference on Architectural Support for Programming Languages and Operating Systems_ , 2019, pp. 1015–1029
  * [21] Siyuan Niu, Adrien Suau, Gabriel Staffelbach and Aida Todri-Sanial  “A Hardware-Aware Heuristic for the Qubit Mapping Problem in the NISQ Era”  In _IEEE Transactions on Quantum Engineering_ 1 IEEE, 2020, pp. 1–14
  * [22] Sutapa Saha, Tamal Guha, Some Sankar Bhattacharya and Manik Banik  “Distributed Computing Model: Classical vs. Quantum vs. Post-Quantum”  In _arXiv preprint arXiv:2012.05781_ , 2020
  * [23] Peter W Shor  “Algorithms for quantum computation: discrete logarithms and factoring”  In _Proceedings 35th annual symposium on foundations of computer science_ , 1994, pp. 124–134  Ieee
  * [24] Yuhong Song et al.  “Dancing along Battery: Enabling Transformer with Run-time Reconfigurability on Mobile Devices”  In _arXiv preprint arXiv:2102.06336_ , 2021
  * [25] Francesco Tacchino et al.  “Quantum implementation of an artificial feed-forward neural network”  In _Quantum Science and Technology_ 5.4 IOP Publishing, 2020, pp. 044010
  * [26] Bochen Tan and Jason Cong  “Optimal layout synthesis for quantum computing”  In _2020 IEEE/ACM International Conference On Computer Aided Design (ICCAD)_ , 2020, pp. 1–9  IEEE
  * [27] Swamit S Tannu and Moinuddin K Qureshi  “A case for variability-aware policies for nisq-era quantum computers”  In _arXiv preprint arXiv:1805.10224_ , 2018
  * [28] Hanrui Wang et al.  “Quantumnas: Noise-adaptive search for robust quantum circuits”  In _arXiv preprint arXiv:2107.10845_ , 2021
  * [29] Zhepeng Wang et al.  “Exploration of Quantum Neural Architecture by Mixing Quantum Neuron Designs”  In _arXiv preprint arXiv:2109.03806_ , 2021
  * [30] Xiaosi Xu et al.  “Variational algorithms for linear algebra”  In _Science Bulletin_ Elsevier, 2021
  * [31] Lei Yang et al.  “Co-exploration of neural architectures and heterogeneous asic accelerator designs targeting multiple tasks”  In _2020 57th ACM/IEEE Design Automation Conference (DAC)_ , 2020, pp. 1–6  IEEE
  * [32] Lei Yang et al.  “Co-exploring neural architecture and network-on-chip design for real-time artificial intelligence”  In _2020 25th Asia and South Pacific Design Automation Conference (ASP-DAC)_ , 2020, pp. 85–90  IEEE