# From ORAN to Cell-Free RAN: Architecture, Performance Analysis, Testbeds and
Trials

Yang Cao, Ziyang Zhang, Xinjiang Xia, Pengzhe Xin, Dongjie Liu, Kang Zheng,
Mengting Lou, Jing Jin, Qixing Wang, Dongming Wang, Yongming Huang, Xiaohu
You, Jiangzhou Wang Y. Cao, Z. Zhang, P. Xin, K. Zheng, D. Wang, Y. Huang and
X. You are with National Mobile Communications Research Laboratory, Southeast
University, Nanjing 210096, China (email: wangdm, xhyu@seu.edu.cn). Y. Cao, X.
Xia, D. Liu, D. Wang, Y. Huang and X. You are also with Purple Mountain Lab.,
M. Lou, J. Jin, Q. Wang are with the China Mobile Research Institute, Beijing
100053, China. J. Wang is with the School of Engineering, University of Kent,
Canterbury CT2 7NT, U.K. (email: j.z.wang@kent.ac.uk).

###### Abstract

Open radio access network (ORAN) provides an open architecture to implement
radio access network (RAN) of the fifth generation (5G) and beyond mobile
communications. As a key technology for the evolution to the sixth generation
(6G) systems, cell-free massive multiple-input multiple-output (CF-mMIMO) can
effectively improve the spectrum efficiency, peak rate and reliability of
wireless communication systems. Starting from scalable implementation of CF-
mMIMO, we study a cell-free RAN (CF-RAN) under the ORAN architecture. Through
theoretical analysis and numerical simulation, we investigate the uplink and
downlink spectral efficiencies of CF-mMIMO with the new architecture. We then
discuss the implementation issues of CF-RAN under ORAN architecture, including
time-frequency synchronization and over-the-air reciprocity calibration, low
layer splitting, deployment of ORAN radio units (O-RU), artificial intelligent
based user associations. Finally, we present some representative experimental
results for the uplink distributed reception and downlink coherent joint
transmission of CF-RAN with commercial off-the-shelf O-RUs.

###### Index Terms:

cell-free massive MIMO, radio access network, open radio access network, 6G

##  I Introduction

With the commercialization of the fifth generation new-radio (5G-NR) networks,
academia and industry have started the research on 5G-Advanced and sixth
generation (6G) technologies [1]. Multiple transmission/reception points
(Multi-TRP) technique which is also called coordinated multiple point (CoMP)
in the forth generation long term evolution (4G LTE) has been considered to be
a key technology for improving spectral efficiency, peak rates and reliability
[2, 3]. Although CoMP was proposed in 4G LTE, it was not widely used in
commercial systems until Release 16 of 5G-NR presented a standardized
implementation of non-coherent Multi-TRP, and currently Release 18 is working
on a standard for coherent joint transmission(CJT) of Multi-TRP [4]. As an
evolution of CoMP, cell-free massive multiple-input multiple-output (CF-mMIMO)
was proposed[5], which also has been considered as a potential 6G technique
[6]. Disregarding the traditional cellular architecture and realizing cell-
free has been the desire of academia and industry for many years. In the study
of CF-mMIMO, the scalable implementation of joint transceiver has the
potential to break the cellular architecture [7].

Network topology and architecture are very essential to the implementation of
cell-free systems. Over the past few years, various kinds of CF-mMIMO
architectures have been presented. Although fully centralized implementation
of CF-mMIMO with a single central processing unit (CPU) has the best
performance, it is not scalable since the fronthaul requirements grow with the
number of access points. Radio stripe architecture proposed by Ericsson [8] is
a sequential approach of distributed processing, and theoretically it is
scalable. In [9], a hybrid architecture was presented, where the scalable CF-
mMIMO network was shaped by several disjoint network-centric clusters (NCC),
and each user was served by several access points belong to different NCC in a
user-centric manner [10]. In [11], a fog massive MIMO architecture was
presented, and multiple edge processing units were introduced to serve their
coordination regions. While most of the current work focuses on the network
topology and transceiver design of CF-mMIMO, how to efficiently implement the
“cell-free” networking from the perspective of the architecture of RAN and its
interfaces requires further study.

Figure 1: Uplink MRC with distributed implementation.

With standardized interfaces, open radio access networks (open RAN or ORAN)
introduce a more open ecology for cellular mobile communications. Currently,
ORAN has developed a rich set of open interfaces for the interoperability of
radio access networks (RANs)[12]. Some operators have deployed ORAN equipments
in commercial 5G networks. The open interfaces enable innovative technologies
to be introduced to improve the performance of cellular mobile communication
systems, while software enabled implementations of baseband processing unit
(BBU) allow for faster development and upgrades [13]. Under ORAN architecture,
some experimental validations of 6G techniques can be developed rapidly with
software implementation. It is also convenient to carry out performance
validation of some key technologies in the existing commercial networks. In
[14] a novel 6G network architecture has been proposed, in which the
hierarchical intelligence, real-time and non-real-time artificial intelligence
(AI) assisted resource allocation have been demonstrated with ORAN-based
testbed. In [15], an edge intelligence-based RAN architecture has been
validated on ORAN 5G networks.

Since ORAN defines a set of open interfaces between radio units (RU) and BBU,
the verification of key technologies of 6G physical layer (PHY) can be
realized by using a generic server with a reference design developed by some
cooperations, such as Intel FlexRAN [13] or Nvdia Aerial [16]. Therefore ORAN
may play an important role in the research and development of CF-mMIMO. On the
one hand, ORAN provides a sophisticated clock and time synchronization
mechanism for multiple RUs. With over-the-air (OTA) reciprocity calibration
[17, 18, 19], downlink CJT has been verified in experimental systems with
commercial off-the-shelf (COTS) ORAN RUs (O-RUs) [20]. On the other hand, the
scalability of CF-mMIMO can be implemented with distributed transceivers,
effectively reducing signalling interactions and significantly reducing the
complexity. The architecture of ORAN provides an excellent platform for the
implementation of cell-free systems. Ranjbar _et al_. [21] put the first
effort to implement CF-mMIMO under the ORAN architecture. Basically, the
architectures in [11, 21] follow the idea of coordination regions, and do not
take full advantage of the favorable propagation effect achieved by CF-mMIMO
system [22], and the scalability is also not discussed. In [23], a general
ORAN based cell-free architecture was further elaborated by introducing the
virtualized RAN functions, and user-centric O-RU selection was also proposed
by using Q-learning algorithm. However, the architectures in [21] and [23]
still follow the current Option 7.2x splitting between O-RU and ORAN based
distributed unit (DU). Nevertheless, the virtualization of DU has been a long-
term challenge since the concept of cloud RAN was proposed. As mentioned in
[21], the inter-DU interface needs to be further specified to achieve higher
spectral efficiency.

In this paper, a cell-free RAN (CF-RAN) is presented under ORAN architecture.
A splitting option between low-PHY and high-PHY is proposed and a new enhanced
common public radio interface (enhanced CPRI or eCPRI) is introduced, so that
the scalability of the cell-free networks can be achieved. With the new
architecture, the performances of spectral efficiency (SE) for both uplink and
downlink of CF-mMIMO systems are studied, and the deployment of O-RUs is
optimized. To achieve user-centric networking, the association between users
and O-RUs is established with the capability of RAN intelligent controller
(RIC). Then the transmission techniques of CF-RAN including OTA reciprocity
calibration, and distributed transceiver are discussed. The CF-RAN testbeds
with COTS O-RUs have been developed and CJTs have been verified.

The paper is organized as follows. Section II introduces a new implementation
of CF-mMIMO system and gives the theoretical analysis of SE. Section III
discusses the design of CF-RAN under ORAN architecture. Section IV presents
the experimental verification results, followed by the conclusions in Section
V.

The notation adopted in this paper conforms to the following convention.
Uppercase and lowercase boldface letters are used to denote matrices and
vectors, respectively. An $M\times M$ identity matrix is denoted by
${\bf{I}}_{M}$. $\left|\cdot\right|$ denotes the absolute value of a scalar.
${\left[\cdot\right]^{\text{T}}}$, ${\left[\cdot\right]^{*}}$ and
${\left[\cdot\right]^{\text{H}}}$ represent the transpose, conjugate and
Hermitian transpose of a vector or a matrix, respectively. ${\rm diag}({\bf
x})$ is a diagonal matrix with $\bf x$ on its diagonal.
$\mathbb{E}\left[\cdot\right]$ represents mathematical expectation. The
distribution of a circularly symmetric complex Gaussian random variable with
zero mean and variance ${\sigma^{2}}$ is denoted as ${\cal C}{\cal
N}\left({0,{\sigma^{2}}}\right)$.

##  II Scalable Implementations of CF-mMIMO systems and SE Analysis

In this section, we will introduce the basic principle of scalable
implementation of a CF-mMIMO system, and present an improved architecture with
a better trade-off between performance and complexity. Under the new
architecture, we analyze the uplink and downlink SE of CF-mMIMO.

###  II-A Basic principle of scalable implementation of CF-mMIMO system

Consider a CF-mMIMO with $L$ O-RUs and $K$ single-antenna user equipments
(UEs). Assuming that each RU is with $N$ antennas, the total number of
antennas in the system is $LN$. We assume that $L$ is large, and $LN\gg K$. At
the $l$th O-RU, the uplink received signal ${{\bf{y}}_{{\rm{UL}},l}}$ can be
expressed as

$\displaystyle{{\bf{y}}_{{\rm{UL}},l}}=\sum\limits_{k=1}^{K}{{{\bf{h}}_{l,k}}\sqrt{{p_{k}}}{s_{k}}}+{{\bf{z}}_{l}}.$
(1)

where ${s_{k}}$ denotes the transmit symbol of the $k$th UE, ${p_{k}}$ denotes
the uplink transmission power of the $k$th UE, ${{\bf{h}}_{l,k}}$ denotes the
$N\times 1$ channel state information (CSI) from the $k$th UE to the $l$th
O-RU, and ${{\bf{z}}_{l}}\sim{{\cal CN}}\left({0,{\sigma_{\rm
UL}^{2}}{{\bf{{I}}}_{N}}}\right)$ denotes additive Gaussian white noise.

Figure 2: Downlink MRT with distributed implementation.

A scalable uplink reception for a CF-mMIMO system with maximum-ratio-combining
(MRC) is shown in Fig.1. Suppose that each O-RU only knows the CSIs from all
UEs to itself. The MRC-based multi-user detection is locally implemented in
each O-RU. At the $l$th O-RU, the detection results of each UE after
quantization can be expressed as

${\hat{s}_{l,k}}={\mathbb{Q}}\left({{\bf{h}}_{l,k}^{\rm
H}{{\bf{y}}_{{\rm{UL}},l}}}\right),$

where ${\mathbb{Q}}\left(\cdot\right)$ denotes a quantization function. Each
O-RU sends the detection results of all UEs to the combining modules to obtain
the final decision of each UE. For downlink transmission, we can use local
maximum-ratio-transmission (MRT) which is shown in Fig.2.

Basically, the scheme is of distributed nature. The advantages of the
distributed implementation are as follows: firstly, distributed coherent
reception and coherent transmission are implemented in each O-RU without any
CSI exchange among O-RUs; secondly, theoretically, even with simple MRC/MRT,
inter-user interference can be eliminated when the number of O-RUs tends to
infinity [5, 3]; thirdly, with fronthaul network, the high-PHY signal
processing of each UE (including combining after detection, coding and
decoding, modulation and demodulation) can be implemented in different
baseband units. Therefore, a virtualized central processing unit (virtualized
CPU or vCPU) can be obtained and its capability expansion is relatively easy;
fourthly, given the number of UEs, the complexity of the system is just
linearly increasing with the number of O-RUs. Then, with deployment of massive
low-cost O-RUs, the inter-user interference can be suppressed with low
complexity.

As seen, the scalable CF-mMIMO is distributed implementation of the joint
transceiver. Basically, the joint transceiver has been split into two physical
modules: coherent receiving/transmitting, and signal combining/distributing.
Theoretically, these two modules can be implemented in a distributed manner,
thus the system can be scaled up infinitely. However, the distributed
implementation of CF-mMIMO has the following problems compared with the
centralized implementation: firstly, joint multi-user detection and precoding
with centralized implementation, such as minimum-mean-squared-error (MMSE)
detection or regularized zero-forcing (RZF) precoding usually have better
performance than distributed implementation, especially for not a large number
of O-RUs; secondly, as shown in Fig.1, each O-RU needs to send the detection
results of all UEs to the next modules, therefore the fronthaul throughput
increases significantly.

###  II-B New Implementation of Scalable CF-mMIMO system

It can be seen that the scalable implementation of CF-mMIMO requires the O-RU
to have the capability of the physical layer processing. For uplink reception,
O-RU has the functions including channel estimation and detection. For
downlink transmission, O-RU should compute multi-user precoding matrix and
perform digital beamforming. This will make RU implementation more difficult.
Currently, O-RUs do not have these capabilities.

In the deployment of existing ORAN, fronthaul multiplexing (FHM) is adopted to
achieve radio frequency (RF) combining of signals from multiple O-RUs or
multiplexing of multiple cells [24]. RF combining has been widely used in 5G
indoor deployment, which can expand coverage but can not achieve multiplexing
gain. Although FHM for multiplexing multiple cells can be used to implement
multi-TRP, it cannot be used in CF-mMIMO to achieve the scalability.

To achieve the scalability of CF-mMIMO and maintain the existing ORAN
deployment architecture, we introduce edge distributed unit (EDU) to replace
FHM [14, 25]. As shown in Fig.3, multiple O-RUs are connected to an EDU, where
the low-PHY can be implemented. For uplink reception, EDU has the functions
including CSI estimation with demodulation reference signal (DMRS) and
sounding reference signal (SRS), joint multi-user or MIMO detection. For
downlink transmission, the EDU performs reciprocity based multi-user
precoding. For uplink transmission, After multi-user detection, the data
streams from different EDUs are combined in vCPU. Since the combination and
the high-PHY processing is user-specific, we call the module as user-centric
distributed unit (UCDU).

Figure 3: A new implementation of Scalable CF-mMIMO system.

Actually, the traditional centralized implementation and fully distributed
implementation are two special cases of the new architecture. When all of the
O-RU are connected to EDU, the new architecture becomes the centralized
implementation. When the number of EDUs is equal to the number of O-RUs, it
becomes fully distributed implementation. Therefore the new architecture can
be considered as a flexible tradeoff between the two implementations.
Theoretically, the proposed architecture is also a special case of the dynamic
cooperative clustering (DCC) in [26]. However, since the O-RUs are fixedly
connected to EDU, the deployment and implementation of the proposed
architecture are easier than fully DCC. Furthermore, as we will show in the
next section, with dynamic UE-O-RU-EDU association, the new architecture can
almost approach the performance of DCC.

In the following, to analyze the SE of CF-mMIMO with the proposed
architecture, we first introduce the system model and channel model.

It is assumed that there are $M$ EDUs and $J$ UCDUs in the system. The channel
vector from the $k$th UE to all O-RUs is denoted as

${{\mathbf{h}}_{k}}={{\left[\mathbf{h}_{k,1}^{\rm T},\mathbf{h}_{k,2}^{\rm
T},\cdots,\mathbf{h}_{k,L}^{\rm T}\right]}^{\rm T}}\in{{\mathbb{C}}^{LN}},$
(2)

which is modeled as

${{\mathbf{h}}_{k}}={{{\mathbf{\Lambda}}^{1/2}_{k}}}{{\mathbf{g}}_{k}},$

where ${{\mathbf{g}}_{k}}$ represents small-scale fading which is modelled by
correlated Rayleigh fading,

${{\mathbf{\Lambda}}_{k}}={\rm
diag}\left({{\lambda}_{k,1}},{{\lambda}_{k,2}},\cdots,{{\lambda}_{k,L}}\right)\otimes{{\mathbf{I}}_{N}},$

indicates large-scale fading, and
${{\lambda}_{k,l}}=d_{k,l}^{-{{\alpha}_{d}}}$, $d_{k,l}$ denotes the distance
between UE $k$ and O-RU $l$, ${{\alpha}_{d}}$ is the path loss exponent. The
correlated Rayleigh fading channel vector ${{\mathbf{g}}_{k}}$ has the
following distribution,
${{\mathbf{g}}_{k}}\sim{{\mathcal{CN}}}\left(\mathbf{0},{{\mathbf{R}}_{k}}\right)$,
in which
${{\mathbf{R}}_{k}}=\operatorname{diag}\left({{\mathbf{R}}_{k,1}},\cdots,{{\mathbf{R}}_{k,L}}\right)\in{{\mathbb{C}}^{LN\times
LN}}$ represents the block diagonal spatial correlation matrix of UE $k$.

###  II-C Uplink SE of CF-mMIMO systems

For the uplink, the signal detection of the $k$th UE can be expressed as:

${{\hat{s}}_{k}}=\sum\limits_{l=1}^{L}{{{{\delta}}_{k,l}}\mathbf{v}_{k,l}^{\rm
H}{{\mathbf{y}}_{\text{UL},l}}},$ (3)

where $\mathbf{v}_{k,l}$ represents a combining vector, ${{\delta}_{k,l}}$ is
the association indicator between the $l$th O-RU and the $k$th UE, that is,
when the $l$th O-RU serves the $k$th UE, ${{\delta}_{k,l}}=1$, otherwise,
${{\delta}_{k,l}}=0$. The association matrix ${\bf D}_{k}$ of the $k$th UE can
be expressed as

${{\mathbf{D}}_{k}}={\rm
diag}\left({{\delta}_{k,1}},{{\delta}_{k,2}},\cdots,{{\delta}_{k,L}}\right)\otimes{{\mathbf{I}}_{N}}.$

Using dynamic association, i.e., designing the matrix ${\bf D}_{k}$ and the
corresponding multiuser detection, we can implement various receivers.

For the fully centralized implementation of CF-mMIMO, CPU has all the
instantaneous CSI of all UEs. Using the use-and-then-forget (UaTF) method in
[27], after joint multi-user detection at CPU, we can express the signal to
interference plus noise ratio (SINR) of the $k$th UE as

$\displaystyle\gamma_{k}^{(\text{UL},\text{c})}=\frac{p_{k}\left|\mathbb{E}\left\\{\mathbf{v}_{k}^{\text{H}}\mathbf{D}_{k}{\mathbf{h}}_{k}\right\\}\right|^{2}}{\sum\limits_{i=1i\neq
k}^{K}p_{i}\mathbb{E}\left\\{|\mathbf{v}_{k}^{\text{H}}\mathbf{D}_{k}{\mathbf{h}}_{i}|^{2}\right\\}+\sigma_{\text{UL}}^{2}\mathbb{E}\left\\{\|\mathbf{D}_{k}\mathbf{v}_{k}\|^{2}\right\\}},$
(4)

where ${{\mathbf{v}}_{k}}={{\left[\mathbf{v}_{k,1}^{\rm
T},\mathbf{v}_{k,2}^{\rm T},\cdots,\mathbf{v}_{k,L}^{\rm T}\right]}^{\rm T}}$.

$\displaystyle\gamma_{k}^{(\text{UL},\text{d})}=\frac{{{p}_{k}}{{\left|\sum\limits_{m=1}^{M}{\mathbb{E}\left({\mathbf{v}}_{{\rm
EDU},k,m}^{\text{H}}{{\mathbf{D}}_{k,m}}{{\mathbf{h}}_{{\rm
EDU},k,m}}\right)}\right|}^{2}}}{\sum\limits_{i=1,i\neq
k}^{K}{{{p}_{i}}}\mathbb{E}{{\left|\sum\limits_{m=1}^{M}{{\mathbf{v}}_{{\rm
EDU},k,m}^{\text{H}}{{\mathbf{D}}_{k,m}}{{\mathbf{h}}_{{\rm
EDU},i,m}}}\right|}^{2}}+\sigma_{\text{UL}}^{2}\mathbb{E}\left(\sum\limits_{m=1}^{M}{\mathbb{E}\|{{\mathbf{D}}_{k,m}}{{\mathbf{v}}_{{\rm
EDU},k,m}}{{\|}^{2}}}\right)}$ (5) $\displaystyle{\bf{v}}_{{\rm
EDU},k,m}={p_{k}}{\left({\sum\limits_{i=1}^{K}{{p_{i}}{{\bf{D}}_{k,m}}\left({{{\hat{\bf{h}}}_{{\rm
EDU},i,m}}\hat{\bf{h}}_{{\rm
EDU},i,m}^{\rm{H}}+{{\bf{C}}_{i,m}}}\right){{\bf{D}}_{k,m}}}+\sigma_{{\rm{UL}}}^{2}{{\bf{I}}}}\right)^{-1}}{{\bf{D}}_{k,m}}{\hat{\bf{h}}_{{\rm
EDU},k,m}}.$ (6)

For the new implementation, after joint multi-user detection at EDU and
combining at UCDU, the SINR of the $k$th UE can be expressed as (5). Takeing
MMSE receiver as an example, the combining vector ${\bf{v}}_{{\rm EDU},k,m}$
at the $m$th EDU is expressed as (6), where the length of combining vector
${\bf{v}}_{{\rm EDU},k,m}$ depends on the total number of antennas in the
$m$th EDU, ${{\hat{\bf{h}}}_{{\rm EDU},i,m}}$ is the channel vector of the
$i$th UE to the $m$th EDU, ${{\tilde{\bf{h}}}_{{\rm EDU},i,m}}$ is the
corresponding channel estimation error vector, and
${\mathbf{C}}_{i,m}=\mathbb{E}\\{{\tilde{\mathbf{h}}_{{\rm
EDU},i,m}\tilde{\mathbf{h}}_{{\rm EDU},i,m}^{\text{H}}}\\}$ is the covariance
matrix of $\tilde{\mathbf{h}}_{{\rm EDU},i,m}$, ${{\mathbf{D}}_{k,m}}$ is the
association matrix. Obviously, the fully distributed implementation is a
special case for $M=L$, and the centralized implementation can be considered
as a special case for $M=1$.

###  II-D Downlink SE of CF-mMIMO systems

For downlink transmission of a time-division duplexing (TDD) system, downlink
CSIs can be obtained by uplink sounding and reciprocity calibration. With CSIs
at transmitter, CJT can be realized in CF-mMIMO system. The received signal of
the $k$th UE is given by

$y_{{\text{DL}},k}=\sum\limits_{l=1}^{L}\mathbf{h}_{k,l}^{\text{H}}\left(\sum\limits_{i=1}^{K}{{{\delta}_{i,l}}{{\mathbf{w}}_{i,l}}{s_{i}}}\right)+{{z}_{k}},$
(7)

where ${s_{i}}\in\mathbb{C}$ is the downlink transmission symbol of the $i$th
UE, ${{\mathbf{w}}_{i,l}}$ is the precoding vector of the $l$th O-RU for the
$i$th UE, $z_{k}\sim{\mathcal{CN}}(0,\sigma_{\rm DL}^{2})$.

Similar to uplink transmission, for fully centralized processing, the
precoding vectors are jointly computed at CPU. Suppose that each UE only uses
statistical CSI to obtain the detection of the transmitted symbol. The SINR of
the $k$th UE is given by[7]

$\displaystyle\gamma_{k}^{(\text{DL},\text{c})}=\text{
}\frac{{{\left|\mathbb{E}\left(\mathbf{h}_{k}^{\text{H}}{{\mathbf{D}}_{k}}{{\mathbf{w}}_{k}}\right)\right|}^{2}}}{\sum\limits_{i=1}^{K}{\mathbb{E}}{{\left|\mathbf{h}_{k}^{\text{H}}{{\mathbf{D}}_{i}}{{\mathbf{w}}_{i}}\right|}^{2}}-{{\left|\mathbb{E}\left(\mathbf{h}_{k}^{\text{H}}{{\mathbf{D}}_{k}}{{\mathbf{w}}_{k}}\right)\right|}^{2}}+\sigma_{\text{DL}}^{2}},$
(8)

where ${{\mathbf{w}}_{k}}={{\left[\mathbf{w}_{k,1}^{\rm
T},\mathbf{w}_{k,2}^{\rm T},\cdots,\mathbf{w}_{k,L}^{\rm T}\right]}^{\rm T}}$
is the joint precoding vector.

$\displaystyle\gamma_{k}^{(\text{DL},\text{d})}=\text{
}\frac{{{\left|\sum\limits_{m=1}^{M}{\mathbb{E}}\left({\mathbf{h}}_{{\rm
EDU},k,m}^{\text{H}}{{\mathbf{D}}_{k,m}}{{\mathbf{w}}_{{\rm
EDU},k,m}}\right)\right|}^{2}}}{\sum\limits_{i=1}^{K}{\mathbb{E}}\left|\sum\limits_{m=1}^{M}{{\mathbf{h}}_{{\rm
EDU},k,m}^{\text{H}}{{\mathbf{D}}_{i,m}}{{\mathbf{w}}_{{\rm
EDU},i,m}}}\right|^{2}-\left|\sum\limits_{m=1}^{M}{\mathbb{E}}\left({\mathbf{h}}_{{\rm
EDU},k,m}^{\text{H}}{{\mathbf{D}}_{k,m}}{{\mathbf{w}}_{{\rm
EDU},k,m}}\right)\right|^{2}+\sigma_{\text{DL}}^{2}}$ (9)
$\displaystyle{\bf{w}}_{{\rm
EDU},k,m}={p_{k}}{\left({\sum\limits_{i=1}^{K}{{p_{i}}{{\bf{D}}_{k,m}}\left({{{\hat{\bf{h}}}_{{\rm
EDU},i,m}}\hat{\bf{h}}_{{\rm
EDU},i,m}^{\rm{H}}+{{\bf{C}}_{i,m}}}\right){{\bf{D}}_{k,m}}}+\sigma_{{\rm{UL}}}^{2}{{\bf{I}}}}\right)^{-1}}{{\bf{D}}_{k,m}}{\hat{\bf{h}}_{{\rm
EDU},k,m}}.$ (10)

For the proposed implementation, precoding vectors are computed at each EDU
with the CSIs of the UEs associated to the EDU. Using the UaTF method in [27]
and MMSE precoding, the SINR of the $k$th UE and the precoding vector
${\bf{w}}_{{\rm EDU},k,m}$ are expressed as (9) and (10) respectively, and
${{p}_{k}}$ is the downlink transmission power allocated to the $k$th UE.
Similar to uplink transmission, the centralized implementation can also be
seen as a special case for $M=1$ and the distributed implementation can be
considered as a special case for $M=L$.

Using random matrix theory and the analysis method in [28], one can obtain the
theoretical expressions of the SINR for both uplink and downlink. In Section
IV, we will give the numerical simulation of the SE performance for CF-mMIMO
with various transceivers.

Figure 4: CF-RAN.

##  III The Design of CF-RAN under ORAN Architecture

In this section, based on the scalable implementation of CF-mMIMO, we will
study the RAN architecture to realize the true cell-free networking. We
propose a CF-RAN architecture, which is illustrated in Fig.4. In CF-RAN, O-RUs
are connected to EDUs, which are connected to UCDU through Ethernet switch.
The O-RUs can be implemented with Option 8 or Option 7-2x [24, 29], and a
conventional CPRI or eCPRI interface may be used between O-RU and EDU. The
function modules of high-PHY, L2 and L3 jointly form vCPU. With ORAN E2
interface, some high layer functions of CF-RAN can be implemented by RAN
intelligent controller (RIC).

In the following, we discuss the key technologies of CF-RAN under ORAN
architecture, including physical layer implementations, O-RU and EDU
deployment, dynamic UE-O-RU-EDU association, and dynamic UE-UCDU association.

###  III-A The synchronization issue of CF-RAN

Phase synchronization between O-RUs is very important for downlink CJT. As
shown in [20], even with common reference clock and timing, the phase drift of
each O-RU is different due to the independent local oscillator (LO), and it
will introduce large performance degradation for downlink CJT. However, the
configuration with time and frequency synchronizations is still beneficial
because it can avoid frequency offset between O-RUs which will greatly reduce
the overhead of frequency offset estimation.

Therefore, to avoid frequency offset and frequent OTA reciprocity calibration
among O-RUs, it is preferable to use SYNCE and 1588 precision time protocol
(PTP) to synchronize all the O-RUs in the system, for example, a configuration
which is called lower layer split C2 (LLS-C2) defined in ORAN can be adopted
[24]. A primary reference time clock/telecom grandmaster (PRTC/T-GM) (acting
as SYNCE and PTP master) is implemented in the fronthaul network to distribute
network timing toward EDU, O-RU and UCDU. With full timing support, all
Ethernet switches in the fronthaul function as telecom boundary clock.

To achieve reciprocity based CJT, OTA reciprocity calibration among O-RUs is
also very important. In [20], a design of calibration reference signal (CARS)
for OTA reciprocity calibration was proposed, which can make full use of the
flexible frame structure of 5G-NR, and can be completely transparent to COTS
O-RU and commercial UEs. The results in [20] also showed that the calibration
coefficient changes with the phase drift of the LO’s phase-locked loop. Since
the O-RUs use independent LOs, the phase of the calibration coefficients
between O-RUs change rapidly within $\pm 30^{\circ}$ with respect to its mean.
Then the CJT algorithms should take into account the phase drift of the
calibration coefficients, otherwise it will cause large performance loss.

###  III-B The interface between EDU and UCDU

In CF-RAN, since the traditional physical layer is implemented separately in
EDU and UCDU, the interface between them is important to the implementation.
As shown in the previous section, the splitting of the two modules can be
placed in the MIMO detection and layer mapping. Referring to the
implementation architecture of ORAN, the splitting Option 7-2x is very similar
to this case. However, unlike the Option 7-2x, the EDU in the CF-RAN
undertakes more physical layer signal processing. Therefore, we need to
further study a new Option 7-2x splitting protocol.

In the following, we assume that O-RU follows Option 8. Figure 5 shows the
downlink functional split for various physical layer channels and transmission
modes. In the downlink, orthogonal frequency division multiplexing (OFDM)
phase compensation [30], inverse fast Fourier transform (inverse FFT/IFFT),
cyclic prefix (CP) addition, resource element (RE) mapping, precoding, and
reciprocity calibration based precoding matrix computation functions reside in
the EDU. The rest of the PHY functions including layer mapping, modulation,
scrambling, rate matching and coding reside in the UCDU. The generation of
calibration reference signal (CARS) [20] for OTA reciprocity calibration is
performed at the UCDU.

Figure 5: Lower layer DL split description. Figure 6: Lower layer UL split
description.

The uplink functional split for various physical layer channels and
transmission modes are illustrated in Figure 6. In the uplink, OFDM phase
compensation (for all channels except physical random access channel (PRACH))
[30], FFT, CP removal, RE-demapping, channel estimation as well as timing
advanced (TA) estimation and interference-plus-noise (IpN) estimation with
demodulation reference signal (DMRS), equalization, physical uplink control
channel (PUCCH) extraction, filtering for PRACH, sounding reference signal
(SRS) channel estimation and measurement, CARS extraction functions reside in
the EDU. The rest of the PHY functions including combining, demodulation,
descrambling, rate dematching and decoding, OTA calibration coefficients
computation reside in the UCDU.

For PUSCH, linear MMSE detection can be adopted as a baseline. Some soft
interference cancellation (SIC) detections, such as MMSE-SIC, expectation
propagation detection [31] can also be considered. The detection results of
each layer of each UE are quantized and compressed and sent to the UCDU
together with the SINR. At the UCDU, using weighted combination according to
the SINR, we obtain the final detection results.

Figure 7: The interface between EDU and UCDU.

For the interface messages between UCDU and EDU, we have designed four message
types, as shown in Figure 7: scheduling messages, uplink and downlink data
messages, measurement messages, reciprocity calibration message. The
scheduling message includes time-frequency resource scheduling information of
different physical channels, as well as multi-user grouping information. The
data messages contain IQ data of different physical channels. Different from
the other uplink physical channels, the PUSCH data message contains the data
streams after equalization/detection. Measurement messages are used to report
channel measurement results to UCDU, such as TA, IpN, SINR, SRS channel
measurement, etc. If PRACH detection is performed in EDU, the detection result
of PRACH can be reported to UCDU through measurement message. The received
CARS is sent to UCDU through a specified data message, and the OTA calibration
coefficients are sent to each EDU through a reciprocity calibration message.
The eCPRI packet for these messages can be designed referring to the fronthaul
specification [24].

Remark 1: In this paper, we present an example of the implementation of CF-
RAN. Actually, the functional split can be very flexible. Firstly, O-RU can be
designed according to Option 7-2x, and the LLS-C3 [24] can also be used to
achieve time-frequency synchronization of O-RUs. However, less functions at
O-RU allow limiting the number of required real time calculations as well as
required memory requirement which means lower O-RU complexity and cost.
Secondly, when the number of O-RUs is large, the random access should be
carefully designed for transmitting synchronization signal block (SSB). The
PRACH processing should be further studied, and it could also be performed at
EDU. Nevertheless, placing most functions at UCDU and less functions at EDU
can keep the interface simple and limits the required associated control
messages. Thirdly, when the number of O-RUs is large, linear equalization can
achieve a better tradeoff between complexity and performance. Maximum-
likelihood detection such as sphere decoding can be used, and the log-
likelihood-ratio compression and combining should also be studied to achieve
better overall performance. Fourthly, the deployments with mmWave RU using
only analog beamforming are also possible with the same interface design.

###  III-C The Deployment of O-RU and EDU

To reduce the cost of a practical deployment, adjacent O-RUs are usually
connected to the same EDU. We call this scheme as clustering deployment of
O-RUs and EDUs. However, in CF-RAN system, when a large number of O-RUs and
multiple EDUs are densely deployed in an area, the clustering deployment even
has worst performance. Considering uplink reception, the joint detection
results of multiple EDUs should be further combined to obtain the final
decision. Intuitively, for clustering deployment, only the UEs on the boundary
covered by the multiple EDUs can achieve the cooperative transmission gains.
That is, if we adopt interleaving deployment of O-RUs and EDUs (which is
illustrated in Figure 4), we can obtain more performance gain from cooperative
transmission. Then, we should study the deployment of O-RU and EDU.

Since the locations of O-RUs are usually known, we present the optimized
connection between O-RUs and EDUs by using genetic algorithm (GA). The
objective is to minimize the sum of the distances between the O-RUs of each
EDU. Since each EDU has the equal processing capability, we make the
constraint that the number of O-RUs connected to each EDU should be as even as
possible. We set the population as an $M$-base sequence $\mathcal{L}$ and the
length of $\mathcal{L}$ is the number of O-RUs $L$. When the gene
$\mathcal{L}_{i}=m$ which is the $i$th element of $\mathcal{L}$, the $i$th
O-RU is connected to the $m$th EDU.

We design the fitness function as

$\displaystyle f(x)=$ (11)
$\displaystyle\frac{1}{\sum\limits_{p\in\mathcal{P}}{\sum\limits_{q\in\mathcal{Q}}{\cdots}}\sum\limits_{u\in\mathcal{U}}{\sum\limits_{v\in\mathcal{V}}{\left({{d}_{p,q}^{2}}+\cdots+{{d}_{p,v}^{2}}+\cdots+{{d}_{u,v}^{2}}\right)^{1/2}}}},$

where ${d}_{p,v}$ is the distance between the $p$th and $v$th O-RUs. We assume
that $\mathcal{P},\mathcal{Q},\cdots,\mathcal{U},\mathcal{V}$ as the O-RU set
of the $p,q,\cdots,u,v$th EDU respectively and $\mathcal{T}$ is the complete
set of O-RUs and the constraints can be described as

$\displaystyle\mathcal{P}\cup\mathcal{Q}\cup\cdots\mathcal{U}\cup\mathcal{V}=\mathcal{T},$
(12a)
$\displaystyle\left(\mathcal{P}\cap\mathcal{Q}\right)\cup\left(\mathcal{P}\cap\mathcal{U}\right)\cup\left(\mathcal{P}\cap\mathcal{V}\right)\cup\cdots\cup\left(\mathcal{U}\cap\mathcal{V}\right)=\varnothing,$
(12b) $\displaystyle\|\mathcal{P}|-|\mathcal{Q}\|\leq
1,\|\mathcal{P}|-|\mathcal{V}\|\leq 1,\cdots,\|\mathcal{U}|-|\mathcal{V}\|\leq
1.$ (12c)

The specific implementation steps are shown in Algorithm 1. In the practical
deployment, when both the coverage area and the number of O-RUs are large, we
can use the above algorithm locally first, and then expand it.

Algorithm 1 GA-based interlacing Algorithm

1:Input: Number of O-RU $L$, Distance matrix between O-RUs $\mathbf{D}$,
Crossover rate $c$, Mutation rate $m$, population size $n_{p}$

2:Pop = makingPopulation($n_{p}$) under constraint of (12);

3:$i=0$;

4:while $i<gen$ do

5: Pop.scores = fitnessFunction(Pop) using (11);

6: newPop = [ ];

7: $j=0$;

8: while $j<numberofchildP$ do

9: [child1,child2] = selection(Pop) with the probability of normalized fitness
value

10: Crossover(child1,child2,$c$);

11: Mutation(child1,child2,$m$);

12: if child meet system constraints (12) then

13: newPop.append(child1)

14: newPop.append(child2)

15: $j=j+2$

16: end if

17: end while

18: newPop.scores = fitnessFunction(newPop) using (11)

19: Pop.append(newPop)

20: SortBasedOnScroes(Pop)

21: Pop = Pop$(0:N-1)$

22: $i=i+1$.

23:end while

24:Output: Best O-RU Group $\mathcal{L}$ = minimumScore(Pop).

###  III-D Dynamic UE-O-RU-EDU Association

In practical systems, the effective number of O-RUs communicating with a UE is
limited due to large-scale fading. Accordingly, the number of EDUs associated
with a UE is also limited. To reduce the complexity of the transceiver in EDU
and reduce the fronthaul transmission, we should dynamically select the
serving O-RUs and corresponding EDUs for each UE according to its location.
With dynamic UE-O-RU-EDU association, we can also achieve DCC and user-centric
CF-mMIMO system.

We propose a Q-learning based dynamic UE-O-RU-EDU association strategy, which
can effectively use the storage complexity to reduce the computational
complexity. We assume that UCDU has the statistical CSIs (large-scale fading)
between UEs and O-RUs. According to the SE analysis method in Section II, when
all the O-RUs serve all UEs, we obtain the rate of each UE. We propose the
following Q-learning algorithm to obtain the UE-O-RU-EDU association.

In the Q-learning algorithm, each EDU serves as an agent and the state of each
EDU at time $t$ is
${{s}_{m,t}}=\left[{{\delta}_{1}},{{\delta}_{2}},\cdots,{{\delta}_{K}}\right]$.
When the $k$th UE is served by EDU $m$, ${{\delta}_{k}}=1$, otherwise,
${{\delta}_{k}}=0$. We use a dynamic $\varepsilon$-greedy strategy.

The reward function is defined as

$r^{(t)}={{\chi}_{t}}\times\frac{{{R}_{\text{sum},t}}}{{{R}_{\text{sum,all}}}-{{R}_{\text{sum},t}}}.$

where the fronthaul constraint and the achievable sum rate are comprehensively
considered as follows: ${\chi_{t}}$ represents whether the fronthaul
constraint is satisfied, if so, it is 1, otherwise, it is 0;
${R_{{\rm{sum}},t}}$ represents the sum rate at time $t$, and
${{R_{{\rm{sum,all}}}}}$ represents the sum rate when all EDUs and UEs are
associated.

The update of the Q-table is according to

$\displaystyle Q_{T}[s^{(t)},{a}^{(t)}]$
$\displaystyle=(1-\alpha)Q_{T}[s^{(t)},{a}^{(t)}]$
$\displaystyle+\alpha\left(r^{(t+1)}+\kappa{\rm
max}_{a}\left\\{Q_{T}\left[s^{(t+1)},{a}^{(t+1)}\right]\right\\}\right),$ (13)

where $\alpha$ denotes learning rate, $\kappa$ denotes discount factor.

The Q-table is continuously updated through the rewards generated by each
action execution, so as to obtain the Q-value for taking a specific action in
a specific state; the agent makes an intelligent decision by observing the
environment, that is, whether to associate EDU with the UE. The agent can gain
experience and adjust the action strategy during the training process, select
the action that can get the maximum reward according to the Q-value, and
realize the association between UEs and EDUs.

Algorithm 2 describes the specific implementation steps. The association
algorithm can be implemented in vCPU. Since the UE-O-RU-EDU association is
near real time (near RT) implementation, it can be realized in near-RT RIC.

Algorithm 2 Q-learning based UE-O-RU-EDU association algorithm

1:Input: Learning rate $\alpha$, discount factor $\kappa$, $\rm{environment}$,
attenuation rate parameter $\varphi$, size of the action
set$\left|\text{action}\right|$, ${{\varepsilon}_{\text{init}}}$, rates when
EDUs and UEs are all associated ${{R}_{\text{sum,all}}}$, $EP$

2:$s^{(t)}\leftarrow\mathbf{0}$, $e\leftarrow 0$, $r\leftarrow 0$,
$\varepsilon\leftarrow 0$;

3:$t_{end}\leftarrow\rm{environment.lastServive}$()

4:while $e<EP$ do

5: $Q_{T}\leftarrow\mathbf{0}$

6: $t\leftarrow 0$

7: $r\leftarrow 0$

8: while $t\leqslant t_{end}$ do

9: if $s^{(t)}==\emptyset$ then

10: $s^{(t)}\leftarrow\rm{environment.getState}(t)$

11: end if

12: Using dynamic $\varepsilon$-greedy strategy

13:
$\varepsilon\left(e\right)\leftarrow{{\varepsilon}_{\text{init}}}{{\left(1-{{\varepsilon}_{\text{init}}}\right)}^{\frac{e}{\varphi\times\left|\text{action}\right|}}}$

14: ${{a}^{t}}\leftarrow\mathbf{P}{\\{\varepsilon\left(e\right)\\}}{\rm
max}_{a}{\\{Q_{T}[s^{(t)},{a}^{(t)}]\\}}+\mathbf{P}{\\{1-\varepsilon\left(e\right)\\}}\rm{Random.action}$

15: ${\rm{environment}}\leftarrow{\rm{DeployService}}(s^{(t)},{{a}^{t}})$

16:
($s^{\prime}$,$r^{(t)}$,${{\chi}_{t}}$,${R}_{\text{sum},t}$)$\leftarrow\rm{environment.current}(t)$

17:
$Q_{T}[s^{(t)},{a}^{(t)}]\leftarrow(1-\alpha)Q_{T}[s^{(t)},{a}^{(t)}]+\alpha(r^{(t+1)}+\kappa{\rm
max}_{a}\\{Q_{T}[s^{(t+1)},{a}^{(t+1)}]\\})$

18: $s\leftarrow s^{\prime}$

19:
$r^{(t)}\leftarrow{{\chi}_{t}}\times\frac{{{R}_{\text{sum},t}}}{{{R}_{\text{sum,all}}}-{{R}_{\text{sum},t}}}$

20: $t\leftarrow t+1$

21: end while

22: $e\leftarrow e+1$

23:end while

24:Output: UE associated O-RUs action ${\\{a_{t}\\}}_{0}^{t_{end}}$.

###  III-E Dynamic load balance in UCDU

In CF-RAN, the processing of a UE’s data is associated with only one UCDU.
Generally, when the UE-O-RU-EDU association is determined, the UE-UCDU
association can be obtained according to the connection between UCDU and EDU.
In addition, since the implementation of UCDU can be cloudized, the
association between UE and UCDU can be performed by using dynamic load
balancing algorithm [32].

##  IV Numerical and Experimental results

In this section, we will give performance evaluation of the CF-RAN with
numerical simulation. We also present the performance of CF-RAN with an
experimental system by using offline evaluations.

###  IV-A Simulation Parameters and Transceiver Schemes

####  IV-A1 Simulation Parameters

It is assumed that the O-RUs are deployed in urban environments with
traditional 2GHz band with the large-scale fading coefficient (channel gain)
in dB as [27]

$\beta_{k,l}(\textrm{dB})=-30.5-36.7\log_{10}\left({d_{k,l}}\right)+F_{k,l},$
(14)

where ${d_{k,l}}$ is the distance (in meter) between O-RU $l$ and UE $k$ and
$F_{k,l}$ is the shadow fading.

The spatial correlation matrix ${\bf R}_{k,l}$ depends on the angular
distribution of the multipath components and we assume the $(m,n)$th element
of ${\bf R}_{k,l}$ can be computed as [28]

$\left[\mathbf{R}\right]_{m,n}=\beta_{k,l}\int\int
e^{j\pi(m-n)\sin(\bar{\varphi})\cos(\bar{\theta})}f(\bar{\varphi},\bar{\theta})d\bar{\varphi}d\bar{\theta}$
(15)

where $\bar{\varphi}$, $\bar{\theta}$ is the azimuth angle and the elevation
angle of a multipath component respectively, $f(\bar{\varphi},\bar{\theta})$
is the joint probability density function of $\bar{\varphi}$ and
$\bar{\theta}$. Specific simulation parameters are shown in the following
table I.

TABLE I: Simulation parameters

[b]  Simulation Parameters Values Uplink transmission power $200\;\rm{mW}$
Antenna height 10 m Area size $200\times 200\;\rm{m}^{2}$ The number of O-RUs,
$L$ 100 Number of antennas per O-RU, $N$ 4 Number of UEs, $K$ 24 Number of
orthogonal pilots,${L_{P}}$ 24 Transmission bandwidth, $B$ 20 MHz Carrier
frequency, $f_{c}$ 2 GHz Azimuth angle, $\bar{\varphi}$ (in radians) 15
Elevation angle, $\bar{\theta}$ (in radians) 15 Noise power spectral density,
$N_{0}$ $-174\;\rm{dBm/Hz}$

####  IV-A2 Transceiver Schemes

In the following subsections, we evaluate the system-level performance with
the centralized processing, fully distributed processing, EDU-based
distributed processing for comparison. The centralized joint MMSE transceiver
has the best performance but the highest implementation complexity. Joint
MRC/MRT can be implemented in a distributed manner without performance loss.
The fully distributed implementation of MMSE transceiver is also called local
MMSE (L-MMSE) transceiver. We also consider the transceivers with DCC, where
the UE-O-RU association is optimized by the Q-learning algorithm. To make a
fair comparison, the proposed CF-RAN with EDU uses the same association as
DCC. In DCC system, only a part of O-RUs serve a UE, and then the joint MMSE
transceiver or fully distributed MMSE transceiver or MRC/MRT are also called
partial MMSE (P-MMSE), LP-MMSE, LP-MRC/MRT, respectively [27, 7]. For CF-RAN,
when all the UEs are associated to an EDU, the distributed transceiver at EDU
with MMSE detection/precoding is named as EDU-MMSE. When only a part of UEs
are associated to an EDU, the distributed transceiver at EDU with MMSE
detection/precoding is called EDU-PMMSE. For the convenience of the reader, we
give the transceiver scheme and its abbreviation in the table II.

TABLE II: Transceiver Schemes

[b]  Transceiver Scheme Abbreviation centralized joint MMSE joint MMSE
centralized joint MRC/MRT joint MRC/MRT fully distributed MMSE L-MMSE joint
MMSE processing with DCC P-MMSE fully distributed MMSE with DCC LP-MMSE fully
distributed MRC/MRT with DCC LP-MRC/MRT EDU-based MMSE without DCC EDU-MMSE
EDU-based MMSE with DCC EDU-PMMSE

Power allocation has a significant impact on the performance of CF-mMIMO
system. For uplink transmission, a fixed transmission power scheme is adopted,
and for downlink transmission, we use the following heuristic power allocation
scheme proposed in [28],

$p_{k}=p_{\max}\frac{\left(\sqrt{\sum_{l\in\mathcal{M}_{k}}\lambda_{k,l}}\right)^{-1}\left(\sqrt{\omega_{k}}\right)^{-1}}{\max_{\ell\in\mathcal{M}_{k}}\sum_{i\in\mathcal{D}_{\ell}}\left(\sqrt{\sum_{l\in\mathcal{M}_{i}}\lambda_{i,l}}\right)^{-1}\sqrt{\omega_{i}}},$
(16)

where $p_{\max}$ is the maximum transmit power per O-RU, $\lambda_{k,l}$ is
the large scale information between UE $k$ and O-RU $l$,
$\omega_{k}=\max_{\ell\in\mathcal{M}_{k}}\mathbb{E}\left\\{\left\|\overline{\mathbf{w}}_{k,\ell}^{\prime}\right\|^{2}\right\\}$
is the part of the normalized precoding vector and $\mathcal{M}_{k}$ is the
corresponding serviced O-RU set.

In the corresponding scheme, power allocation is carried out at the CPU, EDU
and O-RU, respectively. Each O-RU meets the maximum transmission power, and
the total power of all the transmission schemes is the same.

###  IV-B Numerical Results of SE Performance

Figure 8: The CDF of uplink sum SE with interleaving deployment and different
number of EDUs. Figure 9: The CDF of downlink sum SE with interleaving
deployment and different number of EDUs. Figure 10: Comparison of uplink SE
with interleaving deployment and clustering deployments. Figure 11: Comparison
of downlink SE with interleaving deployment and clustering deployments. Figure
12: Test scenario of a CF-RAN prototype system.

Fig. 8 and Fig. 9 demonstrate the cumulative distribution function (CDF) of
the uplink and downlink SE performances for the CF-RAN system. In the
simulations, interleaved deployment between EDU and O-RU is adopted by using
genetic algorithm. In the traditional CF-mMIMO system, the joint MMSE has the
best performance, whereas its implementation is not scalable. The fully
distributed MRC is scalable, however it has large performance loss. The
proposed CF-RAN performs joint detection at the EDU, which can achieve a
flexible tradeoff between centralized processing and distributed processing.
As shown in Fig. 8, given the number of O-RUs, as the number of EDUs
increases, the system performance gradually decreases. When the number of EDUs
is 8, the total number of antennas connected to an EDU is about 50, and the
implementation of EDU is feasible under the current hardware capability. As
can be seen from Fig. 8, its performance can reach $78\%$ of the joint MMSE
receiver. Compared with the fully distributed transceiver, L-MMSE and
L-MRC/MRT, CF-RAN with EDU has significant performance gains. From Fig. 8 and
Fig. 9, it can be seen that the performance with DCC can approach the joint
transmission for both traditional CF-mMIMO and the proposed architecture.

Fig. 10 and Fig. 11 demonstrate the performance of CF-RAN with different O-RU
and EDU deployments. It can be seen that compared with clustering deployment,
the interleaved deployment with GA has better SE performance for both uplink
and downlink transmissions.

###  IV-C Prototype Systems and Performance Evaluation

The CF-RAN experimental system operates at 4.9GHz with 16 O-RUs and 12
prototype UEs, which is shown in Figure 12. Each O-RU is with four antennas,
and each UE is with two antennas. We modify the double-symbol Type 2 DMRS of
5G NR to support 24 ports. In the experimental system, the channel estimation
algorithm of DMRS is Wiener interpolation with uniform power delay profile.
The estimated channels are used for demodulation of the OFDM symbols in the
same slot.

Figure 13 and 14 demonstrate the uplink sum SE of the prototype system with
various reception schemes. We collect the channel estimation results of DMRS,
and analyze the uplink sum SE with offline analysis. Figure 13 shows the CDF
of the sum SE and Figure 14 shows the averaged sum SE. Compared with the
numerical results in Figure 8, we have the similar conclusion. In the
prototype system, compare to the total number uplink data streams, the total
number of antennas is not so large. In this scenario, it can be seen that when
the total number of antennas connected to EDU is larger than the total number
of uplink data streams, the performance gain is most significant.

Figure 15 and 16 demonstrate the downlink sum SE of the prototype system. For
downlink transmission, the CARS configuration proposed in [20] is adopted, and
the total lease-square (TLS) based reciprocity calibration algorithm is used.
Since UE assisted calibration can be considered as perfect CSI feedback, joint
MMSE precoding with UE assisted calibration has the best performance. Unless
specified otherwise self-calibration of the O-RUs is used. In the practical
system, inevitably, there is a time delay when using the calibration
coefficients. Figure 15 shows the SE performance without the delay. It can be
seen that with two EDUs, the performance of EDU-MMSE can approach $80\%$ of
joint MMSE precoding with UE assisted calibration, and is also much better
than L-MMSE. Figure 16 demonstrates the performances with 5ms calibration
delay. It means that the reciprocity calibration based precoding works after 5
ms. Similar to the conclusion in [20], due to LO phase drift, the CSIs change
after 5 ms and it has a significant performance loss for joint preocoding. For
EDU-MMSE, the performance loss is about $10\%$. Even so, EDU-based
implementation is still attractive.

Figure 13: CDF of the uplink sum SE with different reception schemes. Figure
14: Averaged SE of the uplink sum SE with different reception schemes. Figure
15: CDF of the downlink sum SE with different transmitter schemes. Figure 16:
CDF of the downlink sum SE with different transmitter schemes and CSI delay.

##  V Conclusions

In this paper, we provided our considerations for the implementation of “cell-
free” system under ORAN architecture. To begin with, we presented the scalable
implementation of CF-mMIMO system and then introduced EDU to achieve the
tradeoff between joint processing capability and scalability. We derived the
spectral efficiency of the new CF-mMIMO system for both uplink and downlink,
and showed that the traditional fully distributed implementation and fully
centralized implementation were the special cases of the results. We further
elaborated a cell-free RAN under ORAN architecture and then provided the
detailed implementations including time-frequency synchronization for downlink
CJT, low layer splitting for the scalable implementation of distributed
transceiver, O-RU deployment with better system-level performance, UE-O-RU-EDU
association for dynamic user-centric networking. We developed an experimental
system by using commercial ORAN devices to demonstrate the feasibility and
performance gain of CF-RAN with EDU-based implementation.

## References

  * [1] X. You, C.-X. Wang, J. Huang, X. Gao, Z. Zhang, M. Wang, Y. Huang, C. Zhang, Y. Jiang, J. Wang _et al._ , “Towards 6G wireless communication networks: Vision, enabling technologies, and new paradigm shifts,” _Science China Information Sciences_ , vol. 64, no. 1, pp. 1–74, 2021.
  * [2] D. Gesbert, S. Hanly, H. Huang, S. Shamai, O. Simeone, and W. Yu, “Multi-cell MIMO cooperative networks: a new look at interference,” _IEEE Journal on Selected Areas in Communications_ , vol. 28, no. 9, pp. 1–29, 2010.
  * [3] X. You, D. Wang, and J. Wang, _Distributed MIMO and Cell-free Mobile Communication_. Springer, 2021.
  * [4] Samsung, “RP-213015: WID proposal for rel-18 advanced MIMO for NR,” _3GPP TSG-RAN WG1_ , 2021.
  * [5] H. Q. Ngo, A. Ashikhmin, H. Yang, E. G. Larsson, and T. L. Marzetta, “Cell-free massive MIMO versus small cells,” _IEEE Transactions on Wireless Communications_ , vol. 16, no. 3, pp. 1834–1850, 2017.
  * [6] J. Zhang, S. Chen, Y. Lin, J. Zheng, B. Ai, and L. Hanzo, “Cell-free massive MIMO: A new next-generation paradigm,” _IEEE Access_ , vol. 7, pp. 99 878–99 888, 2019.
  * [7] E. Björnson and L. Sanguinetti, “Scalable cell-free massive MIMO systems,” _IEEE Transactions on Communications_ , vol. 68, no. 7, pp. 4247–4261, 2020.
  * [8] G. Interdonato, E. Bj$\ddot{o}$rnson, H. Q. Ngo, P. Frenger, and E. G. Larsson, “Ubiquitous cell-free massive MIMO communications,” _EURASIP J. Wireless Commun. Netw._ , no. 12, p. 197, 2019.
  * [9] G. Interdonato, P. Frenger, and E. G. Larsson, “Scalability aspects of cell-free massive MIMO,” in _IEEE International Conference on Communication (ICC)_ , 2019.
  * [10] C. Pan, M. Elkashlan, J. Wang, J. Yuan, and L. Hanzo, “User-centric c-ran architecture for ultra-dense 5g networks: Challenges and methodologies,” _IEEE Communications Magazine_ , vol. 56, no. 6, pp. 14–20, 2018.
  * [11] A. G. Burr, M. Bashar, and D. Maryopi, “Ultra-dense radio access networks for smart cities: Cloud-RAN, fog-RAN and “cell-free” massive MIMO,” in _International Workshop of CorNer_ , 2018. [Online]. Available: https://arxiv.org/ftp/arxiv/papers/1811/1811.11077.pdf
  * [12] O-RAN Alliance, “O-RAN specifications,” _TS O-RAN_ , 2022. [Online]. Available: https://www.o-ran.org/specifcations
  * [13] Intel Corporation, “Flexrans systems,” _FlexRANs_ , 2022. [Online]. Available: https://www.intel.com/content/www/us/en/developer/topic-technology/edge-5g/tools/flexran.html
  * [14] X. You, Y. Huang, S. Liu, D. Wang _et al._ , “Toward 6G TK$\mu$ extreme connectivity: Architecture, key technologies and experiments,” _Arxiv_ , 2022\. [Online]. Available: https://arxiv.org/abs/2208.01190
  * [15] Y. Liu, Q. Wang, H. Liu, J. Zong, and F. Yang, “Edge intelligence-based RAN architecture for 6G internet of things,” _Discrete Dynamics in Nature and Society_ , 2022. [Online]. Available: https://downloads.hindawi.com/journals/ddns/2022/4955498.pdf
  * [16] Nvdia Corporation, “Nvdia aerial sdk,” _Nvdia Aerial_ , 2022. [Online]. Available: https://developer.nvidia.com/aerial-sdk
  * [17] C. Shepard, H. Yu, N. Anand, E. Li, T. Marzetta, R. Yang, and L. Zhong, “Argos: Practical many-antenna base stations,” in _Proc. 18th Annu. Int. Conf. Mobile Comput. Netw._ , 2012, pp. 53–64.
  * [18] R. Rogalin, O. Y. Bursalioglu, H. Papadopoulos, G. Caire, A. F. Molisch, A. Michaloliakos, V. Balan, and K. Psounis, “Scalable synchronization and reciprocity calibration for distributed multiuser MIMO,” _IEEE Transactions on Wireless Communications_ , vol. 13, no. 4, pp. 1815–1831, 2014\.
  * [19] X. Jiang, A. Decurninge, K. Gopala, F. Kaltenberger, M. Guillaud, D. Slock, and L. Deneire, “A framework for over-the-air reciprocity calibration for TDD massive MIMO systems,” _IEEE Transactions on Wireless Communications_ , vol. 17, no. 9, pp. 5975–5990, 2018.
  * [20] Y. Cao, P. Wang _et al._ , “Experimental performance evaluation of cell-free massive MIMO systems using COTS RRU with OTA reciprocity calibration and phase synchronization,” _Arxiv_ , 2022. [Online]. Available: https://arxiv.org/abs/2208.14048
  * [21] V. Ranjbar, A. Girycki, M. A. Rahman, S. Pollin, M. Moonen, and E. Vinogradov, “Cell-free mmimo support in the o-ran architecture: A phy layer perspective for 5g and beyond networks,” _IEEE Communications Standards Magazine_ , vol. 6, no. 1, pp. 28–34, 2022.
  * [22] Z. Chen and E. Björnson, “Channel hardening and favorable propagation in cell-free massive mimo with stochastic geometry,” _IEEE Transactions on Communications_ , vol. 66, no. 11, pp. 5205–5219, 2018.
  * [23] J. S. Vardakas, K. Ramantas, E. Vinogradov, M. A. Rahman, A. Girycki, S. Pollin, S. Pryor, P. Chanclou, and C. Verikoukis, “Machine learning-based cell-free support in the o-ran architecture: An innovative converged optical-wireless solution toward 6g networks,” _IEEE Wireless Communications_ , vol. 29, no. 5, pp. 20–26, 2022.
  * [24] O-RAN, “Control, user and synchronization plane specification,” _O-RAN Fronthaul Working Group_ , 2020. [Online]. Available: https://www.o-ran.org/specifications
  * [25] D. Wang, X. You, Y. Huang _et al._ , “Full-spectrum cell-free RAN for 6G systems: System design, and experimental results,” _Science China Information Sciences_ , 2023. [Online]. Available: https://www.sciengine.com/SCIS/doi/10.1007/s11432-022-3664-x
  * [26] E. Björnson, N. Jalden, M. Bengtsson, and B. Ottersten, “Optimality properties, distributed strategies, and measurement-based evaluation of coordinated multicell OFDMA transmission,” _IEEE Trans. Signal Process._ , vol. 59, no. 12, pp. 6086–6101, 2011.
  * [27] E. Björnson and L. Sanguinetti, “Making cell-free massive MIMO competitive with MMSE processing and centralized implementation,” _IEEE Transactions on Wireless Communications_ , vol. 19, no. 1, pp. 77–90, 2019.
  * [28] Ö. T. Demir, E. Björnson, L. Sanguinetti _et al._ , “Foundations of user-centric cell-free massive MIMO,” _Foundations and Trends® in Signal Processing_ , vol. 14, no. 3-4, pp. 162–472, 2021.
  * [29] NGMN, “5G RAN CU-DU network architecture, transport options and dimensioning,” _NGMN White Paper_ , 2019. [Online]. Available: https://www.ngmn.org/wp-content/uploads/Publications/2019/190412_NGMN_RANFSX_D2a_v1.0.pdf
  * [30] 3GPP, “Physical channels and modulation,” _3GPP TS 38.211_ , 2019.
  * [31] J. J. Murillo-Fuentes, I. Santos, J. C. Aradillas, and M. Sanchez-Fernandez, “A low-complexity double EP-based detector for iterative detection and decoding in MIMO,” _IEEE Transactions on Communications_ , vol. 69, no. 3, pp. 1538–1547, 2021.
  * [32] K. M. Addali, S. Y. Bani Melhem, Y. Khamayseh, Z. Zhang, and M. Kadoch, “Dynamic mobility load balancing for 5G small-cell networks based on utility functions,” _IEEE Access_ , vol. 7, pp. 126 998–127 011, 2019\.


# A Front Fixing Implicit Finite Difference Method for the American Put
Options Model

Riccardo Fazio111Corresponding author’s e-mail<EMAIL_ADDRESS>Alessandra
Insana and Alessandra Jannelli
Department of Mathematics and Computer Science,
Physical Sciences and Earth Sciences
Viale F. Stagno d’Alcontres 31, I–98166 Messina, Italy

###### Abstract

In this paper, we present an implicit finite difference method for the
numerical solution of the Black-Scholes model of American put options without
dividend payments. We combine the proposed numerical method by using a front
fixing approach where the option price and the early exercise boundary are
computed simultaneously. Consistency and stability properties of the method
are studied. We choose to improve the accuracy of the computed solution via a
mesh refinement based on Richardson’s extrapolation. Comparisons with some
proposed methods for the American options problem are carried out to validate
the obtained numerical results and to show the efficiency of the proposed
numerical methods. Finally, by a posteriori error estimator, we find a
suitable computational grid requiring that the computed solution verifies a
prefixed tolerance.

Key Words: American put options; free boundary problem; front-fixing method;
implicit finite difference scheme; Richardson’s extrapolation.

AMS Subject Classifications: 65M06.

##  1 Introduction

American options are contracts allowing the holder to exercise the option, to
sell or to buy an asset, at a certain price at any time prior to and including
its maturity date. The pricing of American options plays an important role
both in theory and in the real derivative markets. The option pricing model,
developed by Black and Scholes [2] and extended by Merton [17], gives rise to
a partial differential equation governing the value of an option. Schwartz
[24], Brennan and Schwartz [3, 4] were the first to apply a finite difference
method to price American options. The accuracy of their finite difference
method was proved by Jaillet et al. [13]. Other papers present finite
difference methods, (see, for example, Hull and White [12], Duffy [6], Wilmott
et al. [27]). An alternative approach is based on front-tracking methods that
keep track of the boundary and discretize the problem in changing domain; they
have been considered in [10], [19], [25].

When we price an American option we also need to determine the optimal
exercise moment as a function of the value of the underlying asset. This leads
to formulating free boundary of the non-linear problem for the price of the
American option looking for a boundary changing in time to maturity, known as
the optimal exercise boundary. In particular, the American call option problem
is a free boundary problem defined on a finite domain. On the other hand, the
American put option problem is a free boundary problem defined on a semi-
infinite domain so that it is a non-linear model complicated by a boundary
condition at infinity. The difficulty associated with a free boundary problem
can be reduced by using a front-fixing method, an approach relying on a change
of variable to map the changing domain into a fixed domain.

The front-fixing approach has been considered in several papers. In Holmes and
Yang [11], a front-fixing finite element method was used. Tangman et al. [25]
introduced a fourth-order accurate finite difference scheme. In [28], the
original problem is transformed into one more manageable equation with
coefficients not depending on the spatial variable and an explicit update for
the location of the free boundary at each time step was used, and in Zhu et
al. [30], the secant method was employed to solve the non-linear problems.
Moreover, in Tangman et al. [25] and [30], the difference between the prices
of American options and European options was computed. Sevčovič [26] studied
and approximated the early exercise boundary for a class of nonlinear Black-
Scholes equations, applying a fixed-domain transformation and using an
operator splitting iterative numerical technique. Lauko and Ševčovič [16]
introduced a local iterative numerical scheme and compared analytical and
numerical solution for the early exercise boundary position computation of the
American put option. In [29], Zhang and Zhu presented a predictor-corrector
method after the fixed-domain transformation. Nielsen et al. [18] used both
explicit and implicit schemes for solving Black-Scholes model of American
options. While Company et al. [5] presented an explicit finite difference
method for the free boundary value problem under logarithmic front-fixing
transformation. In [7, 8], Fazio proposed one a posteriori error estimator for
the numerical solution for the American put option obtained by an explicit
finite difference scheme.

In this paper, we present an implicit finite difference scheme combined with
the use of the front-fixing method in order to solve the American put option
problem without dividend payments. Preliminary results have been presented to
the International Conference ICNAAM2015 [9]. We use a front-fixing variable
transformation to reformulate the variable domain problem into a non-linear
problem on a fixed rectangular domain. For the non-linear problem, we
introduce a suitable choice of truncated boundary that allows to impose the
asymptotic boundary condition. Then, the original problem is transformed into
a new non-linear partial differential equation where the free boundary appears
as a new variable and computed as part of the solution. We develop an implicit
finite difference method, we investigate the consistency and the stability.
Finally, we choose to validate the obtained numerical results via a mesh
refinement and the Richardson’s extrapolation and we report the comparison
with numerical methods available in the literature.

##  2 The American put options model

Let us suppose that at time $t$ the price of a given underlying asset is $S$.
We consider here the following mathematical model for the value $P=P(S,\tau)$
of an American put option to sell the asset

$\displaystyle{\displaystyle\frac{\partial
P}{\partial\tau}}=\frac{1}{2}\sigma^{2}S^{2}{\displaystyle\frac{\partial^{2}P}{\partial
S^{2}}}+rS{\displaystyle\frac{\partial P}{\partial S}}-rP\quad\mbox{on}\quad
0<\tau\leq T\ ,\ \ S^{*}(\tau)<S<\infty\ ,$ $\displaystyle P(S,0)=\max(E-S,0)\
,\qquad S^{*}(0)=E\ ,$ $\displaystyle\lim_{S\rightarrow\infty}P(S,\tau)=0\ ,$
(1) $\displaystyle P(S^{*}(\tau),\tau)=E-S^{*}(\tau)\
,\qquad{\displaystyle\frac{\partial P}{\partial S}}(S^{*}(\tau),\tau)=-1\ ,$
$\displaystyle P(S,\tau)=E-S\ ,\quad 0\leq S<S^{*}(\tau)\ ,$

where $\tau=T-t$ denotes the time to maturity $T$, $S^{*}(\tau)$ is a free
boundary, that is the unknown early exercise boundary, $\sigma$, $r$ and $E$
are given constant parameters representing the volatility of the underlying
asset, the interest rate and the exercise price of the option, respectively.
Equation (1) is known as the Black-Scholes-Merton equation and it has been
developed by the three economists Fischer Black, Myron Scholes and Robert
Merton in $1973$, [2] and [17].

##  3 The front-fixing method

The front-fixing method is widely employed for solving free boundary problems.
The basic idea of the front fixing method is to use a variable change in order
to remove the free boundary and, then, to transform the original equation into
a new non-linear partial differential equation on a bounded domain, where the
free boundary appears as a new unknown of the problem. The main advantage of
the front-fixing method is that the free boundary is computed directly.

Then, we first transform the Black-Scholes equation into a new parabolic
equation over a bounded domain, we introduce a truncated boundary and, then,
we use finite difference schemes for the new approximate problem over a
bounded domain.

According to Wu and Kwok [28], we consider the dimensionless new variables

$x=\ln\frac{S}{S^{*}(\tau)}\ ,\quad S_{f}(\tau)=\frac{S^{*}(\tau)}{E}\ ,\quad
p(x,\tau)=\frac{P(S,\tau)}{E}\ ,$ (2)

it follows that $S_{f}(\tau)$ is mapped on the fixed line $x=0$, $0\leq
p(x,\tau)\leq 1$ and $0\leq S_{f}(\tau)\leq 1$. Then, the put option problem
(1) can be rewritten as follows

$\displaystyle{\displaystyle\frac{\partial
p}{\partial\tau}}=\frac{1}{2}\sigma^{2}{\displaystyle\frac{\partial^{2}p}{\partial
x^{2}}}+\left(r-\frac{\sigma^{2}}{2}\right){\displaystyle\frac{\partial
p}{\partial
x}}+\frac{1}{S_{f}(\tau)}\frac{dS_{f}}{d\tau}(\tau){\displaystyle\frac{\partial
p}{\partial x}}-r\ p\ ,$ (3) $\displaystyle p(x,0)=0\quad\mbox{for}\quad 0\leq
x\ ,\qquad S_{f}(0)=1\ ,$ (4)
$\displaystyle\lim_{x\rightarrow\infty}p(x,\tau)=0\ ,$ (5) $\displaystyle
p(0,\tau)=1-S_{f}(\tau)\ ,\qquad{\displaystyle\frac{\partial p}{\partial
x}}(0,\tau)=-S_{f}(\tau)\ ,$ (6)

on the domain defined by $0<\tau\leq T$ and $0<x<\infty$.

In order to solve numerically the obtained problem (3-6), we have to consider
a finite computational domain. Then, we introduce a truncated boundary value
$x_{\infty}$, which is a suitable large value where it would be convenient to
impose the asymptotic boundary condition. In other words, we replace the
asymptotic boundary condition (5) with the side condition

$p(x_{\infty},\tau)=0\ .$ (7)

For the choice of $x_{\infty}$ and the accuracy of the related numerical
solution, we can refer to the study by Kangro and Nicolaides [14].

By setting an integer $J$ and a positive value $\mu$, we define the step-sizes

$\Delta x=\frac{x_{\infty}}{J}\ ,\quad\Delta\tau={\mu}{\Delta x}^{2}\ ,$

the integer $N$

$N=\left\lceil\frac{T}{\Delta\tau}\right\rceil\ ,$

where
$\lceil\cdot\rceil:\hbox{I\kern-1.99997pt\hbox{R}}^{+}\rightarrow\hbox{I\kern-1.99997pt\hbox{N}}$
is the ceiling function which maps a real number to the least integer that is
greater than or equal to that number. Therefore, $\mu$ is the grid-ratio

$\mu=\frac{\Delta\tau}{(\Delta x)^{2}}\ .$

Therefore, within the finite domain, we can introduce a mesh of grid-points

$x_{j}=j\Delta x\ ,\quad\tau^{n}=n\Delta\tau\ ,$

for $j=0,1,\dots,J$ and $n=0,1,\dots,N$. We would like to define a numerical
scheme that allows us to compute the grid values

$p_{j}^{n}\approx p(x_{j},\tau^{n})\ ,$

for $j=0,1,\dots,J$, $n=0,1,\dots,N$ and the free boundary values

$S_{f}^{n}\approx S_{f}(\tau^{n})\ ,$

for $n=0,1,\dots,N$.

##  4 A new implicit finite difference scheme

Here we present an implicit finite difference scheme for (3-6).

$\displaystyle\frac{p_{j}^{n+1}-p_{j}^{n}}{\Delta\tau}$
$\displaystyle=\frac{1}{2}\sigma^{2}\frac{p_{j-1}^{n+1}-2p_{j}^{n+1}+p_{j+1}^{n+1}}{(\Delta
x)^{2}}+$ (8)
$\displaystyle+\left(r-\frac{\sigma^{2}}{2}\right)\frac{p_{j+1}^{n+1}-p_{j-1}^{n+1}}{2\Delta
x}+\frac{1}{S_{f}^{n+1}}\frac{S_{f}^{n+1}-S_{f}^{n}}{\Delta\tau}\frac{p_{j+1}^{n+1}-p_{j-1}^{n+1}}{2\Delta
x}-rp_{j}^{n+1}$

for $j=1,2,\dots,J-1$ and $n=0,1,\dots,N-1$.
The initial conditions (4) are

$p_{j}^{0}=0\ ,\quad S_{f}^{0}=1\ ,$

for $j=0,1,\dots,J$. From the boundary conditions (7), we get

$p_{J}^{n+1}=0\ ,$

for $n=0,1,\dots,N$.

Moreover, by considering the governing differential equation (3) at $x_{0}=0$,
$\tau>0$ and taking into account the side conditions (6), one gets a new
boundary condition

$\frac{\sigma^{2}}{2}\frac{\partial^{2}p}{\partial
x^{2}}(0^{+},\tau)+\frac{\sigma^{2}}{2}S_{f}(\tau)-r=0\ ,$ (9)

(see Wu and Kwok [28], Zhang and Zhu [29] or Kwok [15, p. 341]), and its
central finite difference discretization

$\frac{\sigma^{2}}{2}\frac{p_{-1}^{n+1}-2p_{0}^{n}+p_{1}^{n+1}}{(\Delta
x)^{2}}+\frac{\sigma^{2}}{2}S_{f}^{n+1}-r=0\ .$ (10)

From the two boundary conditions (6), using a central finite difference
formula at time $n+1$ and considering (10), we obtain

$p_{0}^{n+1}=1-S_{f}^{n+1},\quad p_{1}^{n+1}=1+\frac{r(\Delta
x)^{2}}{\sigma^{2}}-\left(1+\Delta x+\frac{(\Delta
x)^{2}}{2}\right)S_{f}^{n+1}$ (11)

where $x_{-1}=-\Delta x$ is a fictitious point out of the computational
domain. Considering $\mu={\Delta\tau}/{(\Delta x)^{2}}$ and rearranging the
(8), our implicit numerical scheme can be written as follows

$\bar{a}^{n+1}p_{j-1}^{n+1}+\bar{b}p_{j}^{n+1}+\bar{c}^{n+1}p_{j+1}^{n+1}=p_{j}^{n}$
(12)

for $j=1,2,\dots,J-1$ and $n=0,1,\dots,N-1$, where

$\displaystyle\bar{a}^{n+1}$
$\displaystyle=\frac{\mu}{2}\left[-\sigma^{2}+\left(r-\frac{\sigma^{2}}{2}\right)\Delta
x\right]+\frac{1}{S_{f}^{n+1}}\frac{S_{f}^{n+1}-S_{f}^{n}}{2\Delta x}$
$\displaystyle\bar{b}$ $\displaystyle=1+\mu\sigma^{2}+r\Delta\tau$ (13)
$\displaystyle\bar{c}^{n+1}$
$\displaystyle=\frac{\mu}{2}\left[-\sigma^{2}-\left(r-\frac{\sigma^{2}}{2}\right)\Delta
x\right]-\frac{1}{S_{f}^{n+1}}\frac{S_{f}^{n+1}-S_{f}^{n}}{2\Delta x}\ .$

Considering (11) and putting $j=1$ in (12), we get

$\bar{c}^{n+1}p_{2}^{n+1}=p_{1}^{n}-\bar{a}^{n+1}(1-S_{f}^{n+1})-\bar{b}\left[1+r\frac{(\Delta
x)^{2}}{\sigma^{2}}-\left(1+\Delta x+\frac{1}{2}(\Delta
x)^{2}\right)S_{f}^{n+1}\right]$ (14)

For $j=2$ in (12) we obtain

$\bar{b}p_{2}^{n+1}+\bar{c}^{n+1}p_{3}^{n+1}=p_{2}^{n}-\bar{a}^{n+1}\left[1+r\frac{(\Delta
x)^{2}}{\sigma^{2}}-\left(1+\Delta x+\frac{1}{2}(\Delta
x)^{2}\right)S_{f}^{n+1}\right]\ .$ (15)

Putting $j=J$ in (12)

$\bar{a}^{n+1}p_{J-1}^{n+1}+\bar{b}p_{J}^{n+1}=p_{J}^{n}\ .$ (16)

For $j=3,4,\cdots,J-1$ we have the equations

$\bar{a}^{n+1}p_{j-1}^{n+1}+\bar{b}p_{j}^{n+1}+\bar{c}^{n+1}p_{j+1}^{n+1}=p_{j}^{n}\
.$ (17)

Then, at each time step we obtain a system $J$ equations, given by (14)-(17),
in $J$ unknowns, $p_{2}^{n+1},p_{3}^{n+1},\cdots,p_{J}^{n+1}$ and
$S_{f}^{n+1}$. The system (14)-(17) can be written in the following compact
form

$A(S_{f}^{n+1})p^{n+1}=f(S_{f}^{n+1})\ ,$ (18)

where $p^{n+1}=(p_{2}^{n+1},p_{3}^{n+1},\cdots,p_{J}^{n+1})$, the coefficients
matrix $A=A(S_{f}^{n+1})\in\mathbb{R}^{J,J-1}$ has following form

$\displaystyle
A(S_{f}^{n+1})=\left(\begin{array}[]{cccccc}\bar{c}^{n+1}&&&&&0\\\
\bar{b}&\bar{c}^{n+1}\\\ \bar{a}^{n+1}&\bar{b}&\bar{c}^{n+1}\\\
&\ddots&\ddots&\ddots&\\\ &&\ddots&\ddots&\ddots&\\\
&&&\bar{a}^{n+1}&\bar{b}&\bar{c}^{n+1}\\\
0&&&&\bar{a}^{n+1}&\bar{b}\end{array}\right)\ ,$

and the mapping $f=f(S_{f}^{n+1}):\mathbb{R}\rightarrow\mathbb{R}^{J}$ is

$\displaystyle
f(S_{f}^{n+1})=\left(\begin{array}[]{l}p_{1}^{n}-\bar{a}^{n+1}(1-S_{f}^{n+1})-\bar{b}\left[1+r\frac{(\Delta
x)^{2}}{\sigma^{2}}-\left(1+\Delta x+\frac{1}{2}(\Delta
x)^{2}\right)S_{f}^{n+1}\right]\\\
p_{2}^{n}-\bar{a}^{n+1}\left[1+r\frac{(\Delta
x)^{2}}{\sigma^{2}}-\left(1+\Delta x+\frac{1}{2}(\Delta
x)^{2}\right)S_{f}^{n+1}\right]\\\ p_{3}^{n}\\\ \vdots\\\
p_{J}^{n}\end{array}\right).$

The system (18) can now be written as a non-linear problem in the form

$F(p^{n+1},S_{f}^{n+1})=A(S_{f}^{n+1})p^{n+1}-f(S_{f}^{n+1})=0\ .$ (21)

The implicit discretization leads to a system of non-linear equations (21) for
the price and the location of the free boundary at each time step.

##  5 Consistency and stability

In this section, we discuss the consistency and stability of the implicit
finite difference scheme (8).

###  5.1 Consistency

We write the PDE (3) and the numerical scheme (8) as follows

$\displaystyle L(p,S_{f})=\frac{\partial
p}{\partial\tau}-\frac{1}{2}\sigma^{2}\frac{\partial^{2}p}{\partial
x^{2}}-\left(r-\frac{\sigma^{2}}{2}\right)\frac{\partial p}{\partial
x}-\frac{1}{S_{f}}\frac{\partial S_{f}}{\partial\tau}\frac{\partial
p}{\partial x}+rp=0\ ,$ $\displaystyle
F(p^{n+1}_{j},S^{n+1}_{f})=\frac{p_{j}^{n+1}-p_{j}^{n}}{\Delta\tau}-\frac{1}{2}\sigma^{2}\frac{p_{j-1}^{n+1}-2p_{j}^{n+1}+p_{j+1}^{n+1}}{(\Delta
x)^{2}}+$
$\displaystyle-\left(r-\frac{\sigma^{2}}{2}\right)\frac{p_{j+1}^{n+1}-p_{j-1}^{n+1}}{2\Delta
x}-\frac{1}{S_{f}^{n+1}}\frac{S_{f}^{n+1}-S_{f}^{n}}{\Delta\tau}\frac{p_{j+1}^{n+1}-p_{j-1}^{n+1}}{2\Delta
x}+rp_{j}^{n+1}=0.$

In order to study the consistency, let us take an arbitrary point $(x,\tau)$
in the domain $[0,x_{\infty}]\times[0,T]$, the mesh point
$(x_{j},\tau^{n+1})$, a numerical finite difference method is consistent with
the differential equation if the local truncation error $T_{j}^{n+1}$, defined
by

$T_{j}^{n+1}(\bar{p},\bar{S}_{f})=F(\bar{p}_{j}^{n+1},\bar{S}_{f}^{n+1})\ ,$

satisfies

$T_{j}^{n+1}(\bar{p},\bar{S}_{f})\longrightarrow 0\qquad as\quad\Delta
x\rightarrow 0\ ,\Delta\tau\rightarrow 0\ ,$

where we denote with $\bar{p}_{j}^{n+1}=p(x_{j},\tau^{n+1})$ the exact
solution value of the PDE and with $\bar{S}_{f}^{n+1}=S_{f}(\tau^{n+1})$ the
exact solution of the free boundary. Using Taylor’s expansion, assuming the
existence of the continuous partial derivatives up to order two in time and up
to order four in space, we obtain

$\displaystyle T_{j}^{n+1}(\bar{p},\bar{S}_{f})=-\frac{\Delta\tau}{2}\
\frac{\partial^{2}p}{\partial\tau^{2}}(x_{j},\tau^{n+1})+O(\Delta\tau)^{2}-\frac{\sigma^{2}}{2}\frac{(\Delta
x)^{2}}{12}\ \frac{\partial^{4}p}{\partial x^{4}}(x_{j},\tau^{n+1})+O(\Delta
x)^{4}$ $\displaystyle-\left(r-\frac{\sigma^{2}}{2}\right)\frac{(\Delta
x)^{2}}{3!}\ \frac{\partial^{3}p}{\partial x^{3}}(x_{j},\tau^{n+1})+O(\Delta
x)^{4}$
$\displaystyle-\frac{1}{S_{f}^{n+1}}\frac{dS_{f}}{d\tau}(\tau^{n+1})\left(\frac{(\Delta
x)^{2}}{3!}\ \frac{\partial^{3}p}{\partial x^{3}}(x_{j},\tau^{n+1})+O(\Delta
x)^{4}\right)$ (22)
$\displaystyle+\frac{1}{S_{f}^{n+1}}\left(\frac{\Delta\tau}{2}\frac{d^{2}S_{f}}{d\tau^{2}}(\tau^{n+1})\
+O(\Delta\tau)^{2}\right)\left(\frac{\partial p}{\partial
x}(x_{j},\tau^{n+1})+\frac{(\Delta x)^{2}}{3!}\ \frac{\partial^{3}p}{\partial
x^{3}}(x_{j},\tau^{n+1})+O(\Delta x)^{4}\right)\ .$

Then, the local truncation error is

$T_{j}^{n+1}(\bar{p},\bar{S}_{f})=O(\Delta\tau)+O(\Delta x)^{2}\ .$ (23)

In our case, we have also to consider the boundary conditions (6) and (9)

$\displaystyle L_{1}(p,S_{f})=$ $\displaystyle p(0,\tau)-1+S_{f}(\tau)=0\ ,$
$\displaystyle L_{2}(p,S_{f})=$ $\displaystyle\frac{\partial p}{\partial
x}(0,\tau)+S_{f}(\tau)=0\ ,$ $\displaystyle L_{3}(p,S_{f})=$
$\displaystyle\frac{\sigma^{2}}{2}\frac{\partial^{2}p}{\partial
x^{2}}(0,\tau)+\frac{\sigma^{2}}{2}S_{f}(\tau)-r=0\ .$

The numerical scheme for the boundary conditions is

$\displaystyle F_{1}(p_{0}^{n+1},S_{f}^{n+1})=p_{0}^{n+1}-1+S_{f}^{n+1}=0\ ,$
$\displaystyle
F_{2}(p_{0}^{n+1},S_{f}^{n+1})=\frac{p_{1}^{n+1}-p_{-1}^{n+1}}{2\Delta
x}+S_{f}^{n+1}=0$ $\displaystyle
F_{3}(p_{0}^{n+1},S_{f}^{n+1})=\frac{\sigma^{2}}{2}\frac{p_{-1}^{n+1}-2p_{0}^{n+1}+p_{1}^{n+1}}{(\Delta
x)^{2}}+\frac{\sigma^{2}}{2}S_{f}^{n+1}-r=0\ .$

Using Taylor’s expansion, the local truncation error for the boundary
conditions is of $O(\Delta x)^{2}$

$\displaystyle T_{1}(\bar{p},\bar{S}_{f})=F_{1}(\bar{p},\bar{S}_{f})=0$
$\displaystyle T_{i}(\bar{p},\bar{S}_{f})=F_{i}(\bar{p},\bar{S}_{f})=O(\Delta
x)^{2}\qquad\text{for}\quad i=2,3\ .$

So, assuming the existence of the continuous partial derivatives up to order
two in time and up to order four in space of the solution of the problem
(3-6), the implicit finite difference method defined by (8) is consistent with
the fixed domain model (3) of order $O(\Delta\tau)+O(\Delta x)^{2}$.

###  5.2 Stability

Now, we perform the Von Neumann analysis to investigate the stability of the
implicit method. The Von Neumann analysis is only valid for linear problems
with constant coefficients, it is not possible to apply it to non-linear
problems or to problems with variable coefficients. In order to apply the
method of Von Neumann, it is necessary to linearize the model and freeze the
coefficients, considering the problem locally [1]. Then, the von Neumann
analysis can be applied, and a stability condition can be derived, this
condition will depend on the frozen coefficients involved.

In our context, the non-linear nature of the model (3) is due to the presence
of the term

$\displaystyle\displaystyle{\frac{1}{S_{f}}\frac{dS_{f}}{d\tau}}.$ (24)

with ${S_{f}}$ unknown of the problem. Then, in the stability analysis, we
decide to take it into account and to derive the stability condition in
relation to the value of $S_{f}$ and of its first derivative. Moreover, note
that we approximate the term (24) with

$\displaystyle\frac{1}{S_{f}^{n+1}}\frac{S_{f}^{n+1}-S_{f}^{n}}{\Delta\tau}\
.$ (25)

that assumes nonpositive values because $S_{f}^{n}$ is positive and not
increasing.

Using the Fourier analysis, we set

$p_{j}^{n+1}=\lambda p_{j}^{n}\qquad p_{j\pm 1}^{n+1}=\lambda e^{\pm ik\Delta
x}p_{j}^{n}$

where $\lambda=\lambda(k\Delta x)$ is called the amplification factor.
Substituting these expressions in the numerical method (12) and dividing by
$p_{j}^{n}$, we obtain

$\lambda=\frac{1}{\bar{a}^{n+1}e^{-ik\Delta
x}+\bar{b}+\bar{c}^{n+1}e^{ik\Delta x}}\ .$

By applying the Euler’s formulas and by (4), after some manipulations, we
derive

$\displaystyle\lambda$ $\displaystyle=$
$\displaystyle\frac{1}{1+r\Delta\tau+2\mu\sigma^{2}\sin^{2}{\displaystyle\frac{k\Delta
x}{2}}-i\mu\Delta
x\left[\left(r-{\displaystyle\frac{\sigma^{2}}{2}}\right)+{\displaystyle\frac{1}{S_{f}^{n+1}}}{\displaystyle\frac{S_{f}^{n+1}-S_{f}^{n}}{\Delta\tau}}\right]\sin(k\Delta
x)}\ .$

Now if we compute the magnitude of amplification factor $\lambda$, then we
have

$\displaystyle|\lambda|^{2}=\frac{1}{(1+A)^{2}+B^{2}}\ ,$

where

$\displaystyle A=r\Delta\tau+2\mu\sigma^{2}\sin^{2}{\displaystyle\frac{k\Delta
x}{2}}$ $\displaystyle B=\mu\Delta
x\left[\left(r-{\displaystyle\frac{\sigma^{2}}{2}}\right)+{\displaystyle\frac{1}{S_{f}^{n+1}}}{\displaystyle\frac{S_{f}^{n+1}-S_{f}^{n}}{\Delta\tau}}\right]\sin(k\Delta
x)\ .$

Because $A>0$, $(1+A)^{2}>1$ and $B^{2}>0$, thus $|\lambda|<1$. So, we can
conclude that the implicit finite difference method defined by (8) for the
fixed domain problem (3-6) is unconditionally stable.

In figure (1), we show the amplification factor module for implicit method at
different time steps $\tau^{n}$. In fact, in the stability analysis, we have
decided to evaluate the amplification factor $\lambda$ depending on the value
of (25). We observe that the variation of (25) does not influence the
stability of the implicit method.

Figure 1: Implicit method. Amplification factor module $|\lambda|$ for
different values of $n$ with $N=320$, $\mu=20$ and $x_{\infty}=1$.

In figure (2) we show the amplification factor module for different values of
$\mu$ at final time $\tau=T$. We observe that the unconditionally stable
implicit finite difference scheme produces results that are stable for any
value of grid-ratio $\mu$.

We can compare the result reported in figure (2) with the ones reported in
figure (3), where we show the amplification factor module $\lambda$, for
different values of $\mu$, of the explicit method [5]

$\displaystyle\lambda$ $\displaystyle=$
$\displaystyle{1-r\Delta\tau-2\mu\sigma^{2}\sin^{2}{\displaystyle\frac{k\Delta
x}{2}}+i\mu\Delta
x\left[\left(r-{\displaystyle\frac{\sigma^{2}}{2}}\right)+{\displaystyle\frac{1}{S_{f}^{n}}}{\displaystyle\frac{S_{f}^{n+1}-S_{f}^{n}}{\Delta\tau}}\right]\sin(k\Delta
x)}\ .$

In order to evaluate $\lambda$, we use the same assumptions given for the
implicit method. We observe that the explicit finite difference scheme
produces results that are stable only for values of grid-ratio $\mu<26$.

The main limitation of the explicit method is the restriction on the choice of
$\Delta\tau$ and $\Delta x$ that have to be chosen in relation at the values
of the parameters of the model $\sigma$ and $r$, the volatility of the
underlying asset and the interest rate, respectively. On the other hand, the
new implicit method is unconditionally stable, no restriction on step sizes is
required.

Figure 2: Implicit method. Amplification factor module $|\lambda|$ for
different values of $\mu$ at time $\tau=T$. The implicit method results to be
stable for any value of $\mu$. Figure 3: Explicit method. Amplification factor
module $|\lambda|$ for different values of $\mu$ at time $\tau=T$. The
explicit method becomes unstable for $\mu\geq 26$.

##  6 A posteriori error estimator

In this section, we describe a posteriori error estimator for the American put
options problem based on Richardson’s extrapolation. In this way, we can find
a more suitable computational grid requiring that the numerical solution
obtained by the implicit method verifies a prefixed error tolerance.

For a scalar $U$ of interest, either a value of the solution $p_{j}^{n}$ or a
free boundary value $S_{f}^{n}$, the numerical error $e$ can be defined by

$e=u-U\ ,$ (26)

where $u$ is the exact, usually unknown, value. When the numerical error is
caused prevalently by the discretization error and in the case of smooth
enough solutions the global error can be split into a sum of powers of the
inverse of $N$

$u=U_{N}+C_{0}\left(\frac{1}{N}\right)^{q_{0}}+C_{1}\left(\frac{1}{N}\right)^{q_{1}}+C_{2}\left(\frac{1}{N}\right)^{q_{2}}+\cdots\
,$ (27)

where $C_{0}$, $C_{1}$, $C_{2}$, $\dots$ are coefficients that depend on $u$
and its derivatives, but are independent of $N$, and $q_{0}$, $q_{1}$,
$q_{2}$, $\dots$ are the true orders of the discretization error, see
Schneider and Marchi [23] and the references quoted therein. Each $q_{k}$ is
usually a positive integer with $q_{0}<q_{1}<q_{2}<\cdots$ and all together
constitute an arithmetic progression of ratio $q_{1}-q_{0}$. The value of
$q_{0}$ is called the asymptotic order or the order of accuracy of the method
or of the numerical value $U_{N}$. By replacing into equation (27) $N=N_{g}$
and $N=N_{g+1}$ and subtracting, to the second obtained equation the first
times $(1/s)^{q_{0}}$, $s=N_{g+1}/N_{g}$, we get the first extrapolation
formula

$u\approx U_{g+1}+\frac{U_{g+1}-U_{g}}{s^{q_{0}}-1}\ ,$ (28)

that has a leading order of accuracy equal to $q_{1}$. This type of
extrapolation is due to Richardson [20, 21]. Taking into account equation
(28), we can conclude that the error estimator by a first Richardson’s
extrapolation is given by

$e_{r}=\frac{U_{g+1}-U_{g}}{s^{q_{0}}-1}\ ,$ (29)

where $q_{0}$ is the order of the numerical method used to compute the
numerical solutions. Hence, (29) gives the real value of the numerical
solution error without knowledge of the exact solution. In comparison with
(29), a safer error estimator can be defined by

$e_{s}=U_{g+1}-U_{g}\ .$ (30)

Of course, $q_{0}$ can be estimated with the formula

$q_{0}\approx{\displaystyle\frac{\log(|U_{g}-u|)-\log(|U_{g+1}-u|)}{\log(s)}}\
,$ (31)

where $u$ is again the exact solution (or, if the exact solution is unknown, a
reference solution computed with a suitable large value of $N$), and both $u$
and $U_{g+1}$ are evaluated at the same grid-points of $U_{g}$.

Within the above framework, in order to improve the numerical accuracy by
using only a small number of grid-nodes, we can generalize (28) introducing
the following repeated extrapolation formula

$U_{g+1,k+1}=U_{g+1,k}+\frac{U_{g+1,k}-U_{g,k}}{s^{q_{k}}-1}\ ,$ (32)

where $g\in\\{0,1,2,\dots,G-1\\}$, $k\in\\{0,1,2,\dots,G-1\\}$,
$s=N_{g+1}/N_{g}$ is the grid refinement ratio, and $q_{k}$ is the true order
of the discretization error. The formula (32) is asymptotically exact in the
limit as $N_{0}$ goes to infinity if we use uniform grids. We notice that to
obtain each value of $U_{g+1,k+1}$ we need two computed solutions $U$ in two
adjacent grids, namely $g+1$ and $g$ at the extrapolation level $k$. For any
$g$, the level $k=0$ represents the numerical solution of $U$ without any
extrapolation. We recall that the theoretical orders of accuracy of the
numerical values $U_{g,k}$, with $N=N_{g}$ and $k$ extrapolations, verify the
relation

$q_{k}=q_{0}+k(q_{1}-q_{0})$ (33)

valid for $k\in\\{0,1,2,\dots,G-1\\}$.

##  7 Numerical results

In this section, we show numerical results obtained by using the implicit
finite difference scheme. Comparisons with some numerical methods available in
the literature for the American put options problem are carried out.

The parameters considered for solving the American put option problem (3-6)
are the following

$r=0.1\ ,\quad\sigma=0.2\ ,\quad T=1\ ,\quad E=1\ .$ (34)

In order to choose the truncated boundary value $x_{\infty}$, we computed the
numerical solution $S_{f}^{N}$ for three different values of $x_{\infty}$. As
we can see in table 1, the change of these values does not affect much the
numerical solution; for this reason, we decided to set $x_{\infty}=1$.

$x_{\infty}$ | $N=10$ | $N=20$ | $N=40$
---|---|---|---
$1$ | $0.8710513685210$ | $0.8661003514438$ | $0.86351373597835$
$2$ | $0.8710513685385$ | $0.8661003514444$ | $0.86351373597828$
$4$ | $0.8710513685388$ | $0.8661003514438$ | $0.86351373597789$
Table 1: Free boundary value $S_{f}^{N}$ at $\tau=T$ for different truncated
boundary locations $x_{\infty}$.

In figure 4, we show the plots of $p_{j}^{N}$ versus $x_{j}$ and of
$S_{f}^{n}$ versus $\tau^{n}$; these results are obtained by the implicit
finite difference scheme, considering the parameters (34) with $N=320$ and
$\mu=20$. In figure (5), the $3$D plot of numerical solution $p_{j}^{n}$ is
shown.

Figure 4: Numerical results obtained by implicit scheme (8): on the left
$p_{j}^{N}$ versus $x_{j}$, on the right, $S_{f}^{n}$ versus $\tau^{n}$,
obtained with $J=80$, $N=320$ and $\mu=20$. Figure 5: Numerical results
obtained by implicit scheme (8), $p_{j}^{n}$ with $J=80$, $N=320$ and
$\mu=20$.

In order to solve the non-linear system, we use the MATLAB routine “fsolve”.
Initially, we used the classical Newton iteration method, but the obtained
numerical results showed less accurate when compared with ones obtained by
Matlab routine _fsolve_. By our numerical experiments, we conclude that the
MATLAB routine, principally regarding the choice of the initial iterate,
reveals more robust than the Newton iteration method and then more suited for
our application.

Looking at the results listed in Table 1, we realize also that, fixing a value
of the truncated boundary $x_{\infty}$, the computed values of $S_{f}^{N}$ for
different values of the grid-steps are in agreement only up to the first
decimal place. For this reason, we decide to improve the numerical accuracy by
performing a mesh refinement applying repeated Richardson’s extrapolation. The
implicit difference scheme is accurate of first-order in time and second-order
in space, then, we can perform a mesh refinement keeping constant the grid-
ratio $\mu$, so that we end up with second-order truncation error
$T_{j}^{n}=O(\Delta\tau^{2})$ in time. Therefore, the global error is of
first-order, which is the value of $q_{0}=1$. In our case the sequence of
$s^{q_{k}}$, for $k=0,1,\dots$, with $s=4$ and $q_{k}=k+1$ is given by $4$,
$16$, $64$, $256$, $1024$, $\dots$.

$N$ | $U_{g,0}$ | $U_{g,1}$ | $U_{g,2}$ | $U_{g,3}$ | $U_{g,4}$ | $U_{g,5}$
---|---|---|---|---|---|---
$5$ | $0.884069$ |  |  |  |  |
$20$ | $0.866100$ | $0.860111$ |  |  |  |
$80$ | $0.863100$ | $0.862100$ | $0.862232$ |  |  |
$320$ | $0.862719$ | $0.862592$ | $0.862625$ | $0.862631$ |  |
$1280$ | $0.862717$ | $0.862716$ | $0.862724$ | $0.862726$ | $0.862726$ |
$5120$ | $0.862738$ | $0.862746$ | $0.862748$ | $0.862748$ | $0.862748$ | $0.862748$
Table 2: Implicit method: Richardson’s repeated extrapolations for the free
boundary value $S_{f}^{N}$ at $\tau=T$.

The numerical results, obtained by repeated Richardson’s extrapolations, for
the values $S_{f}^{N}$ computed with the implicit method are reported in Table
2. The last extrapolated is $U_{5,5}=0.862748$, so we can consider this value
as our benchmark $S_{f}^{N}\approx 0.862748$. Our result can be compared with
the value $S_{f}^{N}\approx 0.862762$ found by Fazio [7], with
$S_{f}^{N}\approx 0.86222$ computed by Nielsen et al. [18] and with
$S_{f}^{N}\approx 0.8628$ found by the explicit method of Company et al. [5].
In Table 3, we compare the option price $P(S,T)$ obtained by different methods
with the following parameters

$\displaystyle T=3\qquad\sigma=0.2\qquad r=0.08\qquad E=100\ .$ (35)

We report the “true value” as the reference offered in [22], the penalty
method (PM) of Nielsen et al. given in [22], the explicit method (EM) of
Company et al. [5] with $\Delta x=0.02$, the explicit method with the
Richardson’s extrapolation (EMR) of Fazio [7] with $\mu=5$ and our implicit
method (IM) without extrapolation, setting $\Delta x=0.02$ and $\mu=5$, and
with Richardson’s extrapolation (IMR) with $\mu=5$. In order to find the value
of the option price $P(S,T)$ in correspondence to each different asset price
$S$ we use a piecewise cubic spline interpolation.

Asset | True | PM | EM | EMR | IM | IMR |
---|---|---|---|---|---|---|---
price S | value | Nielsen [18] | Company [5] | Fazio [7] |  |  |
$90$ | $11.6974$ | $11.7207$ | $11.7054$ | $11.7706$ | $11.6926$ | $11.7707$ |
$100$ | $6.9320$ | $6.9573$ | $6.9309$ | $6.9313$ | $6.9243$ | $6.9315$ |
$110$ | $4.1550$ | $4.1760$ | $4.1564$ | $4.1288$ | $4.1467$ | $4.1291$ |
$120$ | $2.5102$ | $2.5259$ | $2.5151$ | $2.5061$ | $2.5028$ | $2.5064$ |
Table 3: Comparison of American put option price $P(S,T)$ with parameters
(35).

Our goal is to solve the American put option problem with a given tolerance
$\epsilon$, where $0<\epsilon\ll 1$; then we use the error estimator defined
by equation (29), or alternatively by equation (30). To this end, we should
solve the given problem twice, for two grids defined with given values of
$J_{g}=J$ and $J_{g+1}=2J$ of space intervals but for the same value of the
grid-ratio $\mu$. The corresponding time intervals $N_{g}$ and $N_{g+1}$
verify the relation $s=N_{g+1}/N_{g}$. Hence, we can apply (component-wise) to
$p^{n}$ and to $S_{f}^{n}$ the error estimator formula (29), or (30). Then, we
can verify whether, for $n=1,2,\dots,N$,

$\|e_{r}(p^{n})\|_{\infty}\leq\epsilon\qquad|e_{r}(S_{f}^{n})|\leq\epsilon\ .$
(36)

If the two inequalities (36) hold true, for $n=1,2,\dots,N$, then we can
accept the numerical solution computed on the grid defined by $J_{g+1}$ and
$N_{g+1}$, otherwise, we have to increase these two integers and repeat the
computation.

Fig. 6 shows the error estimator results computed by setting $\epsilon=0.005$.
We set $\mu=20$ and start with $J_{0}=5$ and $J_{1}=10$ repeating the
computation by doubling the number of spatial grid-intervals if the required
accuracy is not achieved. Our algorithm stops when $J_{5}=160$ that for
$\mu=20$ means $N_{5}=1280$.

Figure 6: Numerical estimated errors $e_{r}(p^{n})$ and $e_{r}(S_{f}^{n})$
versus $\tau^{n}$ for the explicit (left) and implicit method (right), setting
$\epsilon=0.005$, $\mu=20$, $J_{5}=160$ and $N_{5}=1280$.

From the numerical results shown in Figure 6 we can easily realize that the
greatest errors are found within a few time steps for both numerical methods.
The error of the implicit method is almost twice of the error of the explicit
scheme only in the first time steps, but it becomes smaller in the next time
steps, so we can conclude that the numerical results obtained by the implicit
method are more accurate than those of the explicit one. We suggest that in
order to obtain better accuracy also in the first time steps it can be
developed an adaptive version of the both finite difference schemes.

##  8 Concluding Remarks

In this paper, we consider an American put option model, a free boundary
problem defined on an infinite domain. We overcome the numerical difficulty of
solving a free boundary problem using a front fixing approach combined with
finite difference schemes. We propose an unconditionally stable implicit
finite difference scheme and we prove the consistency. In order to improve the
accuracy solution, we use Richardson’s extrapolation. Comparisons with some
methods available in the literature are carried out to validate the obtained
numerical results. Finally, by a posteriori error estimator, we find a
suitable computational grid requiring that the computed solution verifies a
prefixed tolerance.

Acknowledgement. The research of this work was partially supported by the
University of Messina and by the GNCS of INDAM.

## References

  * [1] T. Aslak and W. Ragnar.  Introduction to Partial Differential Equations: A Computational Approach.  Sringer-verlag New York, 1998.
  * [2] F. Black and M. Scholes.  The pricing of options and corporate liabilities.  J. Polit. Econ., 81:637–654, 1973.
  * [3] M. Brennan and E. Schwartz.  The valuation of American put options.  J. Finance, 32:449–462, 1977.
  * [4] M. Brennan and E. Schwartz.  Finite difference methods and jump processes arising in the pricing of contingent claims: a synthesis.  J. Financ. Quant. Anal., 13:461–474, 1978.
  * [5] R. Company, V. N. Egorova, and L. Jódar.  Solving American option pricing models by the front fixing method: Numerical analysis and computing.  Abstr. Appl. Anal., 2014:Article ID 146745, 9 pages, 2014.
  * [6] D. J. Duffy.  Finite Difference Methods in Financial Engineering: A Partial Differential Equation Approach.  John Wiley & Sons, 2006.
  * [7] R. Fazio.  A posteriori error estimator for a front-fixing finite difference scheme for american options.  In S.I. Ao, L. Gelman, D.W.L. Hukins, A. Hunter, and A.M. Korsunsky, editors, Proceedings World Congress on Engineering 2015, held in London, 1-3 July, volume II, pages 691–696. IAENG, Kwun Tong – Hong Kong, 2015.
  * [8] R. Fazio.  American put option: Richardson’s extrapolation and a posteriori error estimator for a front-fixing finite difference scheme.  Accepted for publication by: Appl. Math. E-Notes, 10:–, 2020.
  * [9] R. Fazio, A. Insana, and A. Jannelli.  Front fixing finite difference schemes for American put options model.  In Proceedings of the International Conference of Numerical Analysis and Applied Mathematics (ICNAAM2015), 2016.
  * [10] H. Han and X. Wu.  A fast numerical method for the black-scholes equation of american options.  SIAM J. Numer. Anal., 41:2081–2095, 2003.
  * [11] A. D. Holmes and H. Yang.  A front-fixing finite element method for the valuation of american options.  SIAM J. Sci. Comput., 30:2158–2180, 2008.
  * [12] J. Hull and A. White.  Valuing derivative securities using the explicit finite difference method.  J. Finac. Quant. Anal., 1:87–100, 1990.
  * [13] P. Jaillet, D. Lamberton, and B. Lapeyre.  Variational inequalities and the pricing of American options.  Acta Appl. Math., 21:263–289, 1990.
  * [14] R. Kangro and R. Nicolaides.  Far field boundary conditions for Black-Scholes equations.  SIAM J. Numer. Anal., 38:1357–1368, 2000.
  * [15] Y.-K. Kwok.  Mathematical Models of Financial Derivatives.  Springer, Springer, 2nd edition, 2008.
  * [16] M. Lauko and D. Ševčovič.  Comparison of numerical and analytical approximations of the early exercise boundary of american put options.  The ANZIAM Journal, 51(4):430–448, 2010.
  * [17] R. C. Merton.  The theory of rational option pricing.  Bell J. Econ. Manag. Sci., 4:141–183, 1973.
  * [18] B. F. Nielsen, O. Skavhaug, and A. Tveito.  Penalty and front-fixing methods for the numerical solution of American option problems.  J. Comput. Finance, 5:69–97, 2002.
  * [19] K. N. Pantazopoulos, E. N. Houstis, and S. Kortesis.  Front-tracking finite difference methods for the valuation of american options.  Comput. Econ., 12:255–273, 1998.
  * [20] L. F. Richardson.  The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam.  Proc. R. Soc. London Ser. A, 210:307–357, 1910.
  * [21] L. F. Richardson and J. A. Gaunt.  The deferred approach to the limit.  Proc. R. Soc. London Ser. A, 226:299–349, 1927.
  * [22] F. A. A. E. Saib, Y. D. Tangman, N. Thakoor, and M. Bhuruth.  On some finite difference algorithms for pricing american options and their implementation in mathematica.  In Proceedings of the 11th International Conference on Computational and Mathematical Methods in Science and Engineering (CMMSE ’11), 2011.
  * [23] F. A. Schneider and C. H. Marchi.  On the grid refinement ratio for one-dimensional advection problems with nonuniform grids.  In Proceedings of COBEM 2005, 2005.  18th International Congress of Mechamical Engineering, 8 pages.
  * [24] E. S. Schwartz.  The valuation of warrants: implementing a new approach.  J. Financ. Econ., 4:79–93, 1977.
  * [25] D. Y. Tangman, A. Gopaul, and M. Bhuruth.  A fast high-order finite difference algorithm for pricing american options.  J. Comput. Appl. Math., 222:17–29, 2008.
  * [26] D. Ševčovič.  An iterative algorithm for evaluating approximations to the optimal exercise boundary for a nonlinear black-scholes equation.  Canadian Applied Mathematics Quarterly, 15:77–97, 2007.
  * [27] P. Wilmott, S. Howison, and J. Dewynne.  The Mathematics of Financial Derivatives.  Cambridge University Press, Cambridge, 1995.
  * [28] L. Wu and Y.-K. Kwok.  A front-fixing method for the valuation of American option.  J. Financ. Eng., 6:83–97, 1997.
  * [29] J. Zhang and S. P. Zhu.  A hybrid finite difference method for valuing American puts.  In S.I. Ao, L. Gelman, D.W.L. Hukins, A. Hunter, and A.M. Korsunsky, editors, Proceedings World Congress on Engineering 2009, held in London, 1-3 July, volume II, pages 1131–1136. IAENG, Kwun Tong – Hong Kong, 2009.
  * [30] Y. L. Zhu, B.M. Chen, H. Ren, and H. Xu.  Application of the singularity separating method to american exotic option pricing.  Adv. Comput. Math., 19:147–158, 2003.