\documentclass{tran-l}
%\usepackage[active]{srcltx} % SRC Specials for DVI Searching
% Over-full v-boxes on even pages are due to the \v{c} in author's name
\vfuzz2pt % Don't report over-full v-boxes if over-edge is small

% THEOREM Environments ---------------------------------------------------
 \newtheorem{thm}{Theorem}[subsection]
 \newtheorem{cor}[thm]{Corollary}
 \newtheorem{lem}[thm]{Lemma}
 \newtheorem{prop}[thm]{Proposition}
 \theoremstyle{definition}
 \newtheorem{defn}[thm]{Definition}
 \theoremstyle{remark}
 \newtheorem{rem}[thm]{Remark}
 \numberwithin{equation}{subsection}
% MATH -------------------------------------------------------------------
 \DeclareMathOperator{\RE}{Re}
 \DeclareMathOperator{\IM}{Im}
 \DeclareMathOperator{\ess}{ess}
 \newcommand{\eps}{\varepsilon}
 \newcommand{\To}{\longrightarrow}
 \newcommand{\G}{\mathcal{G}}
 \newcommand{\N}{\mathcal{N}}
 \newcommand{\M}{\mathcal{M}}
 \newcommand{\BOP}{\mathbf{B}}
 \newcommand{\Real}{\mathbb{R}}
 \newcommand{\Field}{\mathcal{F}}
 \newcommand{\RDPlus}{\Real^d_{+}}
 \newcommand{\RPlus}{\Real^d}
 \newcommand{\Prob}{\mathbb{P}}
 \newcommand{\abs}[1]{\left\vert#1\right\vert}
 \newcommand{\set}[1]{\left\{#1\right\}}
 \newcommand{\crs}[1]{\left\langle#1\right\rangle}
 \newcommand{\eqq}{\equiv}
 \newcommand{\essnorm}[1]{\norm{#1}_{\ess}}
 \renewcommand{\baselinestretch}{1.2}
 \usepackage{graphicx}
%%% ----------------------------------------------------------------------
\begin{document}
% Semimartingale Reflected Brownian Motion
\section*{Application: Semimartingale Reflected Brownian Motion}
This section illustrates the application of the algorithm to a class of diffusion processes that
plays a central role in queuing theory: The Semimartingale Reflected Brownian Motions (abbreviated
as SRBM).

These processes have been shown to arise as diffusion limits of open multiclass queueing networks
operating under ``conventional'' Heavy-Traffic conditions: Consider a sequence of queueing networks
indexed by $n$. The key idea of conventional Heavy-Traffic theory is to scale amounts to express
time in multiples of $n$ and space in multiples of $n^{1/2}$ that functional central limit theorem
applies to a properly normalized queue length process while the sequence of traffic intensities
converges to 1 as $n$ increases.

A formal for the SRBM definition is the following:

\begin{defn} Let $S\equiv\RDPlus$ (the positive d-dimensional orthant). Let $\mu$ be a constant vector
in $\RPlus$, $\sigma$ a $d\times d$ non-degenerate covariance matrix (symmetric and strictly
positive definite), and $R$ a $d \times d$ matrix. For each $x\in S$, a SRBM  associated with the
data $(S,\mu,\sigma,R,x)$ is a $\Field_t$-adapted, $d$-dimensional process $Z = (Z(t):t\geq 0)$
defined on some filtered probability space $(\Omega,\Field.\Field_t,\Prob)$ such that:

\begin{itemize}
\item $Z = X + RY$, $\Prob_x$-a.s.,
\item $\Prob_x$-a.s., $Z$ has continuous paths and $Z(t)\in S$ for all $t\geq 0$,
\item $X$ is a d-dimensional Brownian motion with drift vector $\mu$, covariance matrix $\sigma$
and $X(0)=x$. In addition $X(t)-\mu t$ is a $\Field_t$-martingale,
\item $Y$ is an $\Field_t$-adapted d-dimensional process such that under $\Prob$ it satisfies for
each $j=1,\hdots,d:$
\begin{enumerate}
\item $Y(0)=0$
\item $(Y_i(t):t\geq 0)$ is continuous and non-decreasing,
\item $Y_i(t)$ can increase only when $Z$ hits the face $F_i = \{x\in S:x_i=0\}$.
\end{enumerate}
\end{itemize}
\end{defn}

Loosely speaking, SRBM behaves like Brownian motion in the interior of $S$, and is confined to the
orthant by instantaneous ``reflection'' at the boundary faces, where the direction of  reflection
is dictated by the matrix $R$.

The most general condition currently know to ensure existence and uniqueness (in law) of SRBM in
the orthant is that the matrix $R$ is \emph{completely S}.

\begin{defn}
A $d\times d$ matrix R is said to be $\mathcal{S}$ if there exits a d-dimensional vector $u\geq 0$
such that $R u>0$, and to be a completely-$\mathcal{S}$ matrix if each of its principal submatrices
is an $\mathcal{S}$ matrix.
\end{defn}

This completely $\mathcal{S}$ condition is in fact necessary and sufficient. Regarding its
stationary distribution it has been shown that a sufficient condition for its existence and
uniqueness is that all solutions of an associated deterministic Skorohod problem are attracted to
the origin in finite time. Also, $R^{-1}\mu<0$ is known as a necessary condition.

In order to apply the algorithm to this class of processes we need to check that(i) The BAR
condition is sufficient and necessary, and that (ii) we can provide tightness bounds for the
interior and boundary measures.

\subsection*{Bar Condition}
For $f\in C^2_b(S)$, define
\begin{equation}
Lf  \eqq \frac{1}{2} \sum_{i=1}^d\sum_{j=1}^d \sigma_{ij} \frac{\partial^2 f}{\partial x_i\partial
x_j} + \sum_{i=1}^d\mu_i \frac{\partial f}{\partial x_i}
\end{equation}
\begin{equation}
D_i f \eqq R_{\cdot i} \nabla f\, \mbox{  for }x\in F_i\, (i=1,\hdots,d)
\end{equation}

For this process, BAR takes the following form:
\begin{equation}
\int_{\pi}L f d\pi + \sum_{i=1}^d \int_{F_i} D_i f d\nu_i=0 \quad \forall \,f \in
C^2_b(S)\label{BAR}
\end{equation}
where $\pi$ is the stationary distribution for $Z$ associated with boundary measures $\nu_i\,
(i=1,\hdots,d)$.

Necessity of BAR was first derived by Harrison and Williams when the matrix $R$ is Minkowski
($I-R\geq 0$ and $I-R$ is transient), and later by Dai for the \emph{completely S} case.

Sufficiency was proven by Dai-Kurtz through the following theorem:

\begin{thm} Assume  R is a completely-$\mathcal{S}$ matrix. Suppose that $\pi_0$ is a
probability measure on S with support in the interior of S, and $\pi_1,\hdots,\pi_d$ are positive
finite measures with supports on $F_1,\hdots,F_d$ respectively. If they jointly satisfies the basic
adjoint relationship \ref{BAR}, then $\pi_0$ is the stationary distribution for a
$(\sigma,\mu,R)$-SRBM Z.
\end{thm}

\subsection*{Tightness condition}

The use of Lyapunov functions to bound expectations of Markov processes is a widely used technique
(for a clear exposition see Glynn-Zeevi). Given a function $f$ for which one is interested in
computing stationary moments, the key idea is to find a Lyapunov-like function $g$ that satisfies a
certain ``mean'' drift inequality with respect to f, which in turn leads to bounds on the
stationary expectation of the later. This idea has been used to prove the existence of exponential
moments of the interior stationary distribution for the SRBM: Glynn-Zeevi provided an explicit
construction for the case of $R$ being symmetric and positive definite. Budhiraja and Lee have
proven the existence of such function for the \emph{completely S} case.

We prove that same result holds for the boundary measures.
\begin{thm}\label{Thm-Finite-boundary}
Assume R is a completely-$\mathcal{S}$ matrix, and suppose $\pi$ is the stationary distribution for
$Z$ associated with boundary measures $\nu_i\, (i=1,\hdots,d)$ . there exits a vector $\upsilon>0$
such that

\begin{equation*}
\int_{F_i} e^{\upsilon^tx}\, \nu_i(dx) <\infty \quad , i=1,\hdots,d
\end{equation*}
\end{thm}

\begin{proof}

We will need the following proposition proven by Dai.

\begin{prop}\label{prop-bmeasure}
If $\pi$ is the stationary distribution, associated with boundary measures $\nu_i\,
(i=1,\hdots,d)$, then
\begin{equation}
E_{\pi}\left[ \int_0^t f(Z(s)) dY_i(s)\right] = t \int_{F_i}f \nu_i, \quad (i=1,\hdots,d)
\end{equation}
where $\sigma_i$ denote the (d$-$1)-dimensional Lebesgue measure on the face $F_i$.
\end{prop}
Define $\gamma = - R^{-1}\mu > 0$. Now, take a d-dimensional vector $v>0$ such that $R^tv>0$ (the
completely-$\mathcal{S}$ is close under the transpose and inverse operations). Fix $a>0$ and let
$f(x) = \exp(a x^tv)$. Now let $f_n \in C^2_b(S)$ be a positive function that agrees with $f$ for
$x\in B_S(n)\eqq \set{x\in S, \abs{x}\leq n}$. Applying It\^{o}'s formula, we have that
$$f_n(Z(t))-f_n(z(0)) -\int_0^t L f_n(Z(s))ds - \sum_{i=1}^d \int_0^t D_i f_n(Z(s)) dY_i(s)$$
is a $P_x$ martingale. Let $T_n \eqq \inf\set{t\geq 0, \abs{Z(t)}\geq n}$. Since $f_n>0$ we have
$$E_x[f(Z(t\wedge T_n))]\geq E_x[\int_0^{t\wedge T_N} L f(Z(s))ds] + E_x[\sum_{i=1}^d \int_0^{t\wedge T_N} D_i f(Z(s)) dY_i(s)]$$
For this particular function we have $L f(x) = f(x)\left[-a v^t R \cdot \gamma +
\frac{a^2}{2}\sum_{i=1}^d \sum_{j=1}^d \Gamma_{ij} v_i v_j \right]$ and $D_i f(x) = a f(x) [v^t
R_{\cdot i}]$. Taking $a< \frac{2 v^t R \cdot \gamma}{\sum_{i=1}^d \sum_{j=1}^d \Gamma_{ij} v_i
v_j}$ we have

$$E_x[f(Z(t\wedge T_n))]+ \alpha E_x[\int_0^{t\wedge T_N} f(Z(s))ds] \geq  E_x[\sum_{i=1}^d \beta_i \int_0^{t\wedge T_N} f(Z(s)) dY_i(s)]$$
where $\alpha \eqq a v^t R \cdot \gamma - \frac{a^2}{2}\sum_{i=1}^d \sum_{j=1}^d \Gamma_{ij} v_i
v_j >0$ and $\beta_i\eqq v^tR_{\cdot i}>0\quad i=1,\hdots,d$. Since $f>0$ and taking
$\set{f_n}_{n\in \mathbb{N}}$ increasing, the Monotone convergence theorem implies the following:

$$E_x[f(Z(t))]+ \alpha E_x[\int_0^t f(Z(s))ds] \geq  E_x[\sum_{i=1}^d \beta_i \int_0^t f(Z(s)) dY_i(s)]$$

Taking expectations with respect to $\pi$ (and considering $a>0$ small enough for $E_{\pi}[e^{a
v^tx}] $ to exist) we have that

$$E_{\pi}[f]+ \alpha t E_{\pi}[f] \geq  t \sum_{i=1}^d \beta_i \int_{F_i} f \nu_i(ds)$$
where we have used the result on proposition \ref{prop-bmeasure}, the definition of the stationary
distribution and Tonelli's theorem. Dividing by $t$ and taking $t\rightarrow \infty$ we conclude
the following
$$ \int_{F_i} e^{a v^tx} \nu_i(ds) \leq \frac{\alpha}{\beta_i} E_{\pi}[e^{av^tx}]<\infty , \,i=1\hdots,d$$
\end{proof}

These results indicates that moment bounds exists, but we still need practical bounds to be used in
the setting of the algorithm, for tightness and finiteness of the boundary measures.

\subsubsection*{Practical Bounds}
\begin{itemize}

\item \textbf{Finiteness of boundary measures. }we will use a Lyapunov-function argument. Consider a
d-dimensional vector $v>0$ such that $R^tv>0$ and set $f(x) = x^tv$. Proceeding as in proof of
theorem \ref{Thm-Finite-boundary} we have:

$$E_x[f(Z(t))]- E_x[f(x)] + \alpha  =  E_x[\sum_{i=1}^d \beta_i \int_0^t dY_i(s)]$$
where $\alpha \equiv v^t\mu>0$ and $\beta_i\eqq v^tR_{\cdot i}>0\quad i=1,\hdots,d$. Taking
expectations with respect to $\pi$ we have

$$ \int_{F_i} \nu_i(ds) \leq \frac{\alpha}{\beta_i}\quad \,i=1\hdots,d$$

\item \textbf{Tightness of stationary distribution.} Here we will assume that $R$ is symmetric to
facilitate the explicit construction of a simple Lyapunov function. This can be relaxed by means of
a more clever choice of such function. Consider $f(x) = x^tR^{-1}x$. Proceeding as in the proof of
theorem \ref{Thm-Finite-boundary}, we have

$$-f(x)+  E_x[\int_0^t Z(s)^t \gamma ds] = C + E_x[\sum_{i=1}^d \int_0^t Z_i(s) dY_i(s)]$$
where $C = \sum_{i j} \sigma_{ij}R_{ij}$. Notice that $E_x[\sum_{i=1}^d \int_0^t Z_i(s) dY_i(s)] =
0$ by definition of the SRBM. Taking expectations with respect to $\pi$ we have

$$\int_S \langle e, s \rangle \pi(ds) \leq \frac{C}{\min_j\set{\gamma_j}}$$

For the Minkowski case, a bound can be derived directly from lemma 8.4 in Harrison-Williams.

\item \textbf{Tightness of boundary measures.} Take a d-dimensional vector $v>0$ such that $R^tv>0$ and set
the $d \times d$ symmetric matrix $V$ such that  $V_{ij}=v_j v_i>0$. Notice that $V>0$ and $VR>0$.
Let $f(x) = x^tVx$. Proceeding as in the proof of theorem \ref{Thm-Finite-boundary}, we have

$$E_x[f(Z(t))]+  E_x[\int_0^t Z(s)^t V R \gamma ds] \geq  E_x[\sum_{i=1}^d \int_0^t Z(s)^t V R_{\cdot i} dY_i(s)]$$

Taking expectations with respect to $\pi$ we have
$$\int_{F_i} \langle e, s \rangle \nu_i(ds) \leq \frac{\alpha}{\beta_i}\int_S \langle e, s \rangle \pi(ds) ,\, i=1,\hdots,d$$
where $\alpha \eqq \max_j\set{(V R \gamma)_j}>0$ and $\beta_i \eqq \min_j\set{(V R_{\cdot
i})_j}>0,\, i=1,\hdots,d$ We have implicitly used the finiteness of second moments for the
stationary distribution.
\end{itemize}

\subsection*{Algorithm setting}

To apply the algorithm we still need to specify (i) $\{S_n\}$, the sequence of grid for which we
will approximate the BAR condition, and (ii) $H_m$, the sequence of finite dimensional subspaces
approximating $C^2_b(S)$.

\subsubsection*{Grid choice.} Harrison-Williams proved that the stationary distribution for a ``standard'' SRBM has a
separable density function (in the usual Cartesian coordinates) if and only if the covariance
matrix $\sigma$ satisfies the ``skew symmetric'' condition
\begin{equation}
2 \sigma_{ij} = R_{ij}+R_{ji} \quad \mbox{for} \quad i\neq j\label{skew-symmetric}
\end{equation}
In this case, the marginal distribution for coordinate $i$ is exponential with parameter
$2\gamma_i$. More over, the boundary measures are the restriction of the join distribution to the
corresponding faces of $S$. In practice the matrix $R$ does not fulfill this condition, but the
closer it is to fulfill it, the closer the marginal distributions should be to be exponentials.

Recently, Budhiraja and Lee established the finiteness of the moment generating function of the
steady state distribution in a neighborhood of zero, proving the exponential decay of it.

With this in mind we choose the grid to have an ``exponential spacing on the marginal'':
\begin{equation}
S_n = \set{ x\in S | x_i \in\set{[\log(n)-\log(j)]/\lambda_i \quad j=1,\hdots,n} \quad
i=1,\hdots,n}
\end{equation}
For the boundary, we just project this grind on the corresponding face
\begin{equation}
F^n_i = \set{ x\in S_n | x_i =0} \quad i=1,\hdots,n
\end{equation}

\emph{\textbf{Maybe relate this to the choice of a Reference Density in Dai's work.}}

\textbf{Remark.} In practice, the values for $\lambda$ were selected as follows: First we perform a
first run with the algorithm using low values for $\lambda$. Then these values were adjusted
according to the first moments on the marginal distributions.

\subsubsection*{Subspace choice}

For each $m\geq 1$, we choose
\begin{equation}
H_m = \set{f \mbox{polynomial of degree }\leq m}
\end{equation}

The convergence of $H_m\rightarrow C^2_b(S)$ is a well known result (Proper citation).

\subsection*{Numerical Results}

In this section we will compare numerical results from our algorithm with some known analytical
results of particular instances of SRBM. Also, we will compare our algorithm with previous
algorithms using instances for which either simulation or good approximation results are known.

\subsubsection*{A Two-Dimensional SRBM}
Consider a two-dimensional SRBM associated to the following data:
$$ R = \left(\begin{array}{cc}1 & 0\\-1&1\\\end{array}\right) \quad \quad \sigma = \left(\begin{array}{cc}1 & 0\\0&1\\\end{array}\right) \quad \quad \mu = \left(\begin{array}{c}-1 \\0\\\end{array}\right)$$
For this SRBM the stationary condition $\gamma<0$ reduces to $\mu_1<0$. Harrison (CITE) computed a
closed form solution for the stationary distribution density in polar coordinates:
\begin{equation}
p(x,y) = (2|\mu|)^{3/2}/(\pi^{1/2}) r^{-1/2}\exp(\mu_1r(1+\cos(\theta)))\cos(\theta/2)
\end{equation}
where $(x,y) =(r\cos(\theta),r\sin(\theta))$

Without loss of generality consider $\mu_1 =-1$. It can be shown that (CITE Greenberg)
\begin{equation}
E_{\pi}[x] = 0.5 \quad \quad  E_{\pi}[y] = 0.75
\end{equation}
Taking $n=100$ we used our algorithm to get estimations for these moments. The results are shown on
Table \ref{Two-Dim-Table}.
\begin{table}
  \centering
    \begin{tabular}{|c||c|c|c|c|c|c|c|c|} \hline
        $m$ & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 \\\hline \hline
        $E_{\pi}[x]$ & 0.5 & 0.5 & 0.5 & 0.5 & 0.5 & 0.5 & 0.5 & 0.5 \\\hline
        $E_{\pi}[y]$ & 0.793 & 0.821 & 0.764 & 0.768 & 0.762 &0.758 & 0.760 & 0.750 \\\hline
    \end{tabular}

  \caption{Moment estimates Two-dimensional SRBM}\label{Two-Dim-Table}
\end{table}

We see that our algorithm provides good estimations even for low values of $m$. Each one of the
instances shown on Table \ref{Two-Dim-Table} ran in less than $3$ minutes (on a regular Desktop PC)
on a MATLAB implementation of the algorithm. This type of SRBM is the only one (beside the skew
symmetric type) for which the actual distribution is know, so we can use this knowledge to test our
algorithm. Figure \ref{Two-Dim-Figure} shows our marginal distribution estimations for $n=100$ and
$m=6$, compared with the actual marginal distributions (doted line).
\begin{figure}
\includegraphics[width=16cm]{figure1b.eps}
  \caption{Marginal distribution estimations Two-dimensional SRBM}\label{Two-Dim-Figure}
\end{figure}

It is proven (CITE) the one dimensional marginal distribution in $x$ is exponential, however the
same result does not hold for $y$. We see that our algorithm provides good estimations for the
marginal distributions.

\subsubsection*{Symmetric SRBM's}

A standard SRBM is said to be symmetric if its data has the following properties: $\sigma_{ij} =
\rho$ for $1\leq i<j\leq d$, $\mu_i=-1$ for $i=1,\hdots,d$ and $R_{ij}=R_{ji}=-r\leq 0$ for $1\leq
i<j\leq d$. Positiveness of $\sigma$ implies $-1/(d-1)<\rho<1$ and the \emph{completely S}
condition reduces to $r(d-1)<1$. This type of SRBM arises as a Heavy-Traffic limit of a symmetric
generalized Jackson network. Manipulating the BAR condition, Dai (CITE) shows that
\begin{equation}
m_1=\hdots=m_d = \frac{1-(d-2)r+(d-1)r\rho}{2(r+1)}\label{exact-moment-Symmetric}
\end{equation}
where $m_i = E_{\pi}[x_i]$ $i=1,\hdots,d$. We use our algorithm to compute estimates for these
moments for the case $d=2$. Imitating Dai's work, we let $\rho$ range through $\set{-0.9,
-0.5,0.0,0.5,0.9}$ and $r$ range through $\set{0.2, 0.4,0.6,0.8,0.9, 0.95}$. Table
\ref{Symmetric-Table} shows the relative errors between our estimations and the exact values for
$m=6$ and $n=100$. We see that in each case the relative error is lower that $1\%$.

\begin{table}
  \centering
    \begin{tabular}{|c||c|c|c|c|c|} \hline
        $r/\rho$ & -0.9 & -0.5 & 0.0 & 0.5 & 0.9  \\\hline \hline
        $0.2 $ & 1.98e-4 & 1.65e-4 & 1.26e-3 & 2.32e-3 & 6.66e-4  \\\hline
        $0.4 $ & 3.90e-4 & 4.88e-5 & 2.52e-5 & 2.92e-3 & 2.00e-4  \\\hline
        $0.6 $ & 1.76e-5 & 5.97e-4 & 2.05e-4 & 1.06e-6 & 5.70e-4  \\\hline
        $0.8 $ & 1.53e-5 & 4.27e-5 & 2.22e-5 & 1.17e-3 & 2.16e-3  \\\hline
        $0.9 $ & 1.65e-5 & 4.26e-5 & 5.43e-4 & 1.00e-3 & 3.24e-3  \\\hline
        $0.95$ & 4.57e-6 & 3.42e-5 & 2.46e-4 & 5.98e-4 & 4.82e-3  \\\hline
    \end{tabular}

  \caption{Relative errors for Symmetric SRBM}\label{Symmetric-Table}
\end{table}


\subsubsection*{Skew-symmetric SRBM}
Harrison and Williams proved that a standard SRBM has a product form stationary distribution  if
and only if $\gamma<0$ and condition \ref{skew-symmetric} holds. In this, case we know that the
marginal distribution for $x_i$ is exponential with mean $1/(2\gamma_i)$. Consider a
two-dimensional SRBM associated to the following data:
$$ R = \left(\begin{array}{cc}1 & -0.6\\-0.25&1\\\end{array}\right) \quad \quad \sigma = \left(\begin{array}{cc}1 & -0.425\\-0.425&1\\\end{array}\right) \quad \quad \mu = \left(\begin{array}{c}-0.85 \\0\\\end{array}\right)$$
One can check that condition \ref{skew-symmetric} holds, and that $\gamma^t =(-1,-0.25)$. This
implies that $E_{\pi}[x_1] =0.5$ and  $E_{\pi}[x_2] =2$. Taking $n=100$ we used our algorithm to
get estimates for these moments. The results are shown on Table \ref{Skew-Table}.
\begin{table}
  \centering
    \begin{tabular}{|c||c|c|c|c|c|c|c|c|} \hline
        $m$ & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 \\\hline \hline
        $E_{\pi}[x_1]$ & 0.4942 & 0.4943 & 0.4968 & 0.5014 & 0.5005 & 0.5011 & 0.4996 & 0.5002 \\\hline
        $E_{\pi}[x_2]$ & 2.0095 & 2.0094 & 2.0053 & 1.9976 & 1.9990 &1.9981 & 2.0005 & 1.9996 \\\hline
    \end{tabular}

  \caption{Moment estimates Skew-symmetric SRBM}\label{Skew-Table}
\end{table}

We see that our algorithm provides good estimations even for low values of $m$. Each one of the
instances shown on Table \ref{Two-Dim-Table} ran in less than $3$ minutes (on a regular desktop PC)
on a MATLAB implementation of the algorithm. Since we know that each marginal distribution is
exponential, we can compare our estimation marginal distribution against the actual marginal
distribution for each coordinate. Figure \ref{Skew-Figure} shows our marginal distribution
estimations for $n=100$ and $m=6$, compared with the actual marginal distributions (doted line).
\begin{figure}
\includegraphics[width=16cm]{figure1b.eps}
  \caption{Marginal distribution estimations Skew-symmetric SRBM}\label{Skew-Figure}
\end{figure}

We see that our algorithm provides good estimations for the marginal distributions.

\subsubsection*{Suresh and Whitt's Experiments}

This section follows closely the analysis presented in Chapter 4 of Dai's dissertation (CITE).
Consider a network of d queues in tandem. Let $rho_i$ denote the mean service time at station $i$,
$C^2_{s_i}$ denote the squared coefficient of variation of the service time distribution at station
$i$, and $C^2_a$ denote the squared coefficient of variation for the interarrival time
distribution. Using the method proposed by Harrison and Nguyen, Dai shows that the d-dimensional
current workload process can be approximated by a d-dimensional SBRM associated with the data
$(\sigma,R,\mu)$. When $d=2$is shown that

$$ R = \left(\begin{array}{cc}1 & 0\\-\rho_2/\rho_1&1\\\end{array}\right) \quad \quad \sigma = \left(\begin{array}{cc}\rho^2_1(C^2_a+C^2_{s_1}) & -rho_1\rho_2C^2_{s_1}\\\rho^2_2(C^2_{s_1}+C^2_{s_2})& -rho_1\rho_2C^2_{s_1}\\\end{array}\right) \quad \quad \mu = \left(\begin{array}{c}\rho_1-1 \\\rho_2/\rho_1-1\\\end{array}\right)$$
The arrival rate is assumed to be 1, so that $\rho_i$ represents the traffic intensity at station
$i$.

Suresh and Whitt studied a system of two queues in tandem with the purpose of finding the best
order for the queues to minimize the average total waiting time. They consider various variability
parameter triples $C^2_a,C^2_{s_1},C^2_{s_2}$ for all combinations of the traffic intensities
$\rho_1$ and $\rho_2$ in a representative range. For $C^2_{s_1}\neq C^2_{s_2}$ they considered five
variability triples, namely, (0.5,0.5,2.0), (1.0,0.5,8.0), (1.0,2.0,4.0), (4.5,0.5,1.0) and
,(4.0,1.0,4.0). We will refer them as Case 1 to Case 5, respectively. For $C^2_{s_1}= C^2_{s_2}$
they consider two variability triples, namely, (1.0,0.5,0.5) and (1.0,4.0,4.0). We will refer them
as Case 6 and Case 7 respectively. When  $C^2_{s_1}\neq C^2_{s_2}$, for each queue they consider
four values of $\rho_i:$ 0.3, 0.6,0.8,0.9. When  $C^2_{s_1}\neq C^2_{s_2}$, for each queue they
consider five values of $\rho_i:$ 0.1,0.2,0.3,0.6,0.9.

Extensive simulation experiments were conducted in order to obtain estimates for the expected
steady-state waiting times. When $C^2 = 0.5$, the $E_2$ distribution is used. When $C^2 = 1.0$, the
exponential distribution is used. When $C^2 > 1.0$, the $H_2$ distribution with balanced means is
used.

They also compare their simulation results with QNA approximations. Dai use this set of experiments
to test the QNET method. There, simulation and approximation results for the expected waiting time
at the second station are compared.

As expected, we will use these results to compare our algorithm against the QNA and QNET methods.
Table X-Y (appendix) give simulation estimates, QNET estimates, QNA estimates, and our ``LP''
estimates. Table \ref{Suresh-Whitt-Table} summarizes all balanced heavy traffic cases and gives an
overall comparison of QNA estimates, QNET estimates and LP estimates under heavy traffic. There,
quantities shown represent the minimum between absolute relative difference and absolute
difference, with respect to simulation estimates.

\begin{table}
  \centering
    \begin{tabular}{||c|c|c||c|c|c||c|c|c||} \hline
      \multicolumn{3}{||c||}{} & \multicolumn{3}{|c||}{Order I} & \multicolumn{3}{|c||}{Order II}\\
      \hline
      Case & $\rho_1$ & $\rho_2$ & QNET-D & QNA-D & LP-D & QNET-D & QNA-D & LP-D\\ \hline
      1 & 0.9 & 0.9 & 0.01 & 0.01 & 0.00 & 0.02 & 0.13 & 0.01\\
        & 0.8 & 0.8 & 0.08 & 0.08 & 0.08 & 0.06 & 0.08 & 0.07\\\hline
      2 & 0.9 & 0.9 & 0.03 & 0.06 & 0.03 & 0.01 & 0.21 & 0.01\\
        & 0.8 & 0.8 & 0.04 & 0.01 & 0.02 & 0.08 & 0.07 & 0.07\\ \hline
      3 & 0.9 & 0.9 & 0.04 & 0.02 & 0.04 & 0.04 & 0.19 & 0.05\\
        & 0.8 & 0.8 & 0.00 & 0.04 & 0.01 & 0.02 & 0.02 & 0.01\\ \hline
      4 & 0.9 & 0.9 & 0.13 & 0.30 & 0.14 & 0.19 & 0.29 & 0.19\\
        & 0.8 & 0.8 & 0.14 & 0.08 & 0.11 & 0.09 & 0.00 & 0.07\\ \hline
      5 & 0.9 & 0.9 & 0.01 & 0.12 & 0.01 & 0.01 & 0.01 & 0.00\\
        & 0.8 & 0.8 & 0.09 & 0.14 & 0.08 & 0.05 & 0.05 & 0.05\\ \hline
      6 & 0.9 & 0.9 & 0.01 & 0.08 & 0.02 & 0.01 & 0.10 & 0.01\\
        & 0.8 & 0.8 & 0.01 & 0.04 & 0.01 & 0.01 & 0.04 & 0.01\\ \hline
      7 & 0.9 & 0.9 & 0.04 & 0.18 & 0.05 & 0.06 & 0.07 & 0.04\\
        & 0.8 & 0.8 & 0.02 & 0.05 & 0.03 & 0.01 & 0.06 & 0.01\\ \hline
      \multicolumn{3}{||c||}{Average} & 0.05 & 0.09 & 0.04 & 0.05 & 0.09 & 0.04\\ \hline
    \end{tabular}
  \caption{Overall Comparisons with QNA and QNET Approximations in Heavy Traffic}\label{Suresh-Whitt-Table}
\end{table}

We used $n=100$ and $m=6$ for all experiments. Each run took less than 3 minutes on regular Desktop
PC. As expected, QNET and LP methods give ``better'' estimates than the ones from QNA method, under
balanced heavy traffic conditions. Also, LP method performs as good as the QNET method, most of the
time better.

%\input{Appendix}
\end{document}
