**Fractal Dimension Estimation Methods for Biomedical Images**

Antonio Napolitano, Sara Ungania and Vittorio Cannata

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/48760

## **1. Introduction**

160 MATLAB – A Fundamental Tool for Scientific Computing and Engineering Applications – Volume 3

Peninsula, Portugal. *Resources, Conservation and Recycling*, vol. 56, pp. 7-21

countermeasures. *Journal of Environmental Management*, vol. 88, pp. 479–495

*Process*, RWS Publications, Pittsburgh, PA

*Industrial Engineering*, vol. 61, pp. 226-237

*Energy Reviews*, vol. 14, pp. 1569-1579

*Fuzzy Sets and Systems*, vol. 11, pp. 229-241

*Sets and Systems*, vol. 112, pp. 359-404

pp. 584-599

3071

1432

real option valuation model. *Omega*, vol. 35, pp. 247-57

Pires, Ana.; Chang, Ni-Bin. & Martinho, Graça. (2011). An AHP-based fuzzy interval TOPSIS assessment for sustainable expansion of the solid waste management system in Setúbal

Promentilla, Michael Angelo B.; Furuichi, T. ; Ishii, K. & Tanikawa, N. (2008). A fuzzy analytic network process for multi-criteria evaluation of contaminated site remedial

Saaty, T. L. (1996). *Decision-making with Dependence and Feedback: The Analytic Network* 

Shakhsi-Niaei, M.; Torabi, S. A. & Iranmanesh, S. H. (2011). A comprehensive framework for project selection problem under uncertainty and real-world constraints. *Computers &* 

Smith-Perera, Aida.; García-Melón, Mónica.; Poveda-Bautista, Rocío. & Pastor-Ferrando, Juan-Pascual. (2010). A Project Strategic Index proposal for portfolio selection in electrical company based on the Analytic Network Process. *Renewable and Sustainable* 

Srdjevic, Bojan. (2005). Combining different prioritization methods in the-analytic hierarchy

Van Laarhoven, P. J. M. & Pedrycz, W. (1983). A fuzzy extension of Saaty's priority theory.

Wang, J. & Hwang, W. L. (2007). A fuzzy set approach for R&D portfolio selection using a

Wang, J.; Xu,Y. J. & Li, Z. (2009). Research on project selection system of pre-evaluation of engineering design project bidding. *International Journal of Project Management*, vol. 27,

Wang, Y. M.; Elhag, T. M. S. & Hua, Z. S. (2006). A modified fuzzy logarithmic least squares method for fuzzy analytic hierarchy process. *Fuzzy Sets and Systems*, vol. 157, pp. 3055-

Xu, R. (2000). Fuzzy least-squares priority method in the analytic hierarchy process. *Fuzzy* 

Yu, Jing-Rung. & Cheng, Sheu-Ji. (2007). An integrated approach for deriving priorities in analytic network process. *European Journal of Operational Research*, vol. 180, pp. 1427-

Zadeh, L. A. (1965). Fuzzy sets. *Information and Control*, vol. 8, pp. 338-353

process synthesis. *Computers & Operations Research*, vol. 32, pp. 1897-1919

The use of medical images has its main aim in the detection of potential abnormalities. This goal is accurately achieved with the synergy between the ability in recognizing unique image patterns and finding the relationship between them and possible diagnoses. One of the methods used to aid this process is the extrapolation of important features from the images called texture; texture is an important source of visual information and is a key component in image analysis.

The current evolution of both texture analysis algorithms and computer technology made boosted development of new algorithms to quantify the textural properties of an image and for medical imaging in recent years. Promising results have shown the ability of texture analysis methods to extract diagnostically meaningful information from medical images that were obtained with various imaging modalities such as positron emission tomography (PET) and magnetic resonance imaging (MRI). Among the texture analysis techniques, fractal geometry has become a tool in medical image analysis. In fact, the concept of fractal dimension can be used in a large number of applications, such as shape analysis[1] and image segmentation[2]. Interestingly, even though the fact that self-similarity can hardly be verified in biological objects imaged with a finite resolution, certain similarities at different spatial scales are quite evident. Precisely, the fractal dimension offers the ability to describe and to characterize the complexity of the images or more precisely of their texture composition.

## **2. Fractals**

## **2.1. Fractal geometry**

A fractal is a geometrical object characterized by two fundamental properties: *Self-similarity* and *Hausdorff Besicovich dimension*. A self-similar object is exactly or approximately similar to a part of itself and that can be continuously subdivided in parts each of which is (at least approximately) a reduced-scale copy of the whole. Furthermore, a fractal generally

©2012 Napolitano et al., licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. ©2012 Napolitano et al., licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 2 Will-be-set-by-IN-TECH 162 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 3 Fractal Dimension Estimation Methods for Biomedical Images <sup>3</sup>

shows irregular shapes that cannot be simply described by Euclidian dimension, but, fractal dimension (*f d*) has to be introduced to extend the concept of dimension to these objects. However, unlike topological dimensions the fd can take non-integer values, meaning that the way a fractal set fills its space is qualitatively and quantitatively different from how an ordinary geometrical set does.

Nature presents a large variety of fractal forms, including trees, rocks, mountains, clouds, biological structures, water courses, coast lines, galaxies[3]. Moreover, it is possible to construct mathematical objects which satisfy the condition of self-similarity and that present fd (Figure 1).

**Figure 1.** Sierpinski triangle: starting with a simple initial configuration of units or with a geometrical object then the simple seed configuration is repeatedly added to itself in such way that the seed configuration is regarded as a unit and in the new structure these units are arranged with respect to each other according to the same symmetry as the original units in the seed configuration. And so on.

The objects in Figure 1 are self-similar since a part of the object is similar to the whole and the fractal dimension can be calculated by the equation:

$$D = \frac{\log N}{\log S} \tag{1}$$

Hausdorff formulation[3] is based on the construction of a particular measure, *H<sup>D</sup>*

Intuitively we can sum up the construction as follows: let be *A* a fractal and *C*(*r*, *A*) = {*B*1....*Bk*} a complete coverage of *A* consisting of spheres of diameter smaller than a given

*r*→0

of length *r*. The shorter the ruler, the longer the length measured, a paradox known as the

Hence, when *r* → 0 the effective length of *A* is well approximated. Limit for small *r* calculated

*<sup>δ</sup>* <sup>→</sup> 0 and <sup>H</sup>*<sup>D</sup>*

{*inf* ∑ *i δD*

*δ* :

representing the uniform density of the fractal object.

We define the Hausdorff measure as the function *H<sup>D</sup>*

with *ω<sup>D</sup>* volume of the unit sphere in *R<sup>D</sup>* for integer *D*.

for other values of *D*, however, lead to a degenerate *H<sup>D</sup>*

that provides a method to estimate the dimension *DH*. In the uni-dimensional case *D* = 1 we can easily obtain:

*LD*

*H<sup>D</sup>*

H*D*

*<sup>δ</sup>* the *D*-dimensional Hausdorff measure given by Eq. 3. The *course-grained volume* defined by Eq. 3 normally presents a scaling like:

Therefore, *DH* can be defined as the transition point for the function *H<sup>D</sup>*

*DH*(*A*) = *infD*><sup>0</sup> {*H<sup>D</sup>*

*H<sup>D</sup>*

*<sup>δ</sup>* <sup>∼</sup> *<sup>r</sup>*(1−*DH* ) with *<sup>L</sup><sup>D</sup>*

length. Similar arguments can be applied to *<sup>D</sup>* <sup>=</sup> 2; for *<sup>r</sup>* <sup>→</sup> 0 the measure of *<sup>H</sup><sup>D</sup>*

Although the definition of Hausdorff dimension is particularly useful to operatively define the fd, that presents difficulties when implementing it. In fact, determining the lower bound value of all coverings, as defined in Eq. 5, can be quite complex. For example, let's consider the uni-dimensional case, in which we want to compute the fd of a coastline (Koch Curve). According to Eq. 3 in the case of *D* = 1 the coastline length is measured by a ruler of length *r*. Accuracy of the measure increase with decreasing *r*. For *r* → 0 the coastline will have infinite

*<sup>δ</sup>* (*A*) = *ω<sup>D</sup>* lim

We obtain an approximate measurement of *A*, the so-called *course-grained volume*[4].

r that approximate *A*, so *δ<sup>i</sup>* = *δi*(*Bi*) < *r*.

In the one-dimensional case (*D* = 1), *H<sup>D</sup>*

*coastline paradox*[3].

decreasing with *D*:

from which we derive *DH*.

**3. Methods**

with *H<sup>D</sup>*

covering spheres for *A* with *δ* < *r*:

*δ* ,

*<sup>δ</sup>* that identifies the smallest of all the

Fractal Dimension Estimation Methods for Biomedical Images 163

*<sup>δ</sup>* → ∞ (4)

*<sup>δ</sup>* (*A*) = 0} (5)

*<sup>δ</sup>* = *measured length* (7)

*<sup>δ</sup>* <sup>∼</sup> *<sup>r</sup>*(*D*−*DH* ) (6)

*<sup>δ</sup>* monotonically

*<sup>δ</sup>* → 0.

*<sup>δ</sup>* supplies the length of set *A* measured with a ruler

*<sup>i</sup>* } (3)

where *N* is the number of the auto-similar parts in which an object can be subdivided and *S* is the scaling, that is, the factor needed to observe *N* auto-similar parts. According to the Eq.1, the following values are obtained for the Koch fractal and the Sierpinski triangle:

$$D\_{\rm Koch} = \frac{\log 4}{\log 3} \approx 1.26 \qquad \qquad D\_{\rm Sierpinski} = \frac{\log 3}{\log 3} \approx 1.58 \tag{2}$$

In mathematics, no universal definition of fd exists and the several definitions of fd may lead to different results for the same object. Among the wide variety of fd definitions that have been introduced, the Hausdorff dimension *DH* is surely the most important and the most widely used[4]. Such definition can be theoretically applied to every fractal set but has the disadvantage it cannot always be easily determined by computational methods.

## **2.2. Hausdorff dimension** *DH*

Hausdorff dimension *DH* was introduced in 1918 by mathematical Felix Hausdorff [3]. Since many of the technical developments used to compute the Hausdorff dimension for highly irregular sets were obtained by Abram Samoilovitch Besicovitch, *DH* is sometimes called Hausdorff-Besicovitch dimension.

Hausdorff formulation[3] is based on the construction of a particular measure, *H<sup>D</sup> δ* , representing the uniform density of the fractal object.

Intuitively we can sum up the construction as follows: let be *A* a fractal and *C*(*r*, *A*) = {*B*1....*Bk*} a complete coverage of *A* consisting of spheres of diameter smaller than a given r that approximate *A*, so *δ<sup>i</sup>* = *δi*(*Bi*) < *r*.

We define the Hausdorff measure as the function *H<sup>D</sup> <sup>δ</sup>* that identifies the smallest of all the covering spheres for *A* with *δ* < *r*:

$$H\_\delta^D(A) = \omega\_D \lim\_{r \to 0} \{ \inf \sum\_i \delta\_i^D \} \tag{3}$$

with *ω<sup>D</sup>* volume of the unit sphere in *R<sup>D</sup>* for integer *D*.

We obtain an approximate measurement of *A*, the so-called *course-grained volume*[4].

In the one-dimensional case (*D* = 1), *H<sup>D</sup> <sup>δ</sup>* supplies the length of set *A* measured with a ruler of length *r*. The shorter the ruler, the longer the length measured, a paradox known as the *coastline paradox*[3].

Hence, when *r* → 0 the effective length of *A* is well approximated. Limit for small *r* calculated for other values of *D*, however, lead to a degenerate *H<sup>D</sup> δ* :

$$\mathcal{H}\_{\delta}^{D} \to 0 \quad \text{and} \quad \mathcal{H}\_{\delta}^{D} \to \infty \tag{4}$$

Therefore, *DH* can be defined as the transition point for the function *H<sup>D</sup> <sup>δ</sup>* monotonically decreasing with *D*:

$$D\_H(A) = \inf\_{D \ge 0} \left\{ H\_\delta^D(A) = 0 \right\} \tag{5}$$

with *H<sup>D</sup> <sup>δ</sup>* the *D*-dimensional Hausdorff measure given by Eq. 3.

The *course-grained volume* defined by Eq. 3 normally presents a scaling like:

$$H\_\delta^D \sim r^{(D-D\_H)} \tag{6}$$

that provides a method to estimate the dimension *DH*.

In the uni-dimensional case *D* = 1 we can easily obtain:

$$L\_\delta^D \sim r^{(1-D\_H)} \qquad\qquad\text{with}\quad L\_\delta^D = measured\text{ length}\tag{7}$$

from which we derive *DH*.

#### **3. Methods**

2 Will-be-set-by-IN-TECH

shows irregular shapes that cannot be simply described by Euclidian dimension, but, fractal dimension (*f d*) has to be introduced to extend the concept of dimension to these objects. However, unlike topological dimensions the fd can take non-integer values, meaning that the way a fractal set fills its space is qualitatively and quantitatively different from how an

Nature presents a large variety of fractal forms, including trees, rocks, mountains, clouds, biological structures, water courses, coast lines, galaxies[3]. Moreover, it is possible to construct mathematical objects which satisfy the condition of self-similarity and that present

**Figure 1.** Sierpinski triangle: starting with a simple initial configuration of units or with a geometrical object then the simple seed configuration is repeatedly added to itself in such way that the seed

configuration is regarded as a unit and in the new structure these units are arranged with respect to each other according to the same symmetry as the original units in the seed configuration. And so on.

The objects in Figure 1 are self-similar since a part of the object is similar to the whole and the

*<sup>D</sup>* <sup>=</sup> log *<sup>N</sup>*

where *N* is the number of the auto-similar parts in which an object can be subdivided and *S* is the scaling, that is, the factor needed to observe *N* auto-similar parts. According to the Eq.1,

log 3 <sup>≈</sup> 1.26 *DSierpinski* <sup>=</sup> log 3

In mathematics, no universal definition of fd exists and the several definitions of fd may lead to different results for the same object. Among the wide variety of fd definitions that have been introduced, the Hausdorff dimension *DH* is surely the most important and the most widely used[4]. Such definition can be theoretically applied to every fractal set but has the

Hausdorff dimension *DH* was introduced in 1918 by mathematical Felix Hausdorff [3]. Since many of the technical developments used to compute the Hausdorff dimension for highly irregular sets were obtained by Abram Samoilovitch Besicovitch, *DH* is sometimes called

the following values are obtained for the Koch fractal and the Sierpinski triangle:

disadvantage it cannot always be easily determined by computational methods.

log *<sup>S</sup>* (1)

log 3 <sup>≈</sup> 1.58 (2)

fractal dimension can be calculated by the equation:

*DKoch* <sup>=</sup> log 4

**2.2. Hausdorff dimension** *DH*

Hausdorff-Besicovitch dimension.

ordinary geometrical set does.

fd (Figure 1).

Although the definition of Hausdorff dimension is particularly useful to operatively define the fd, that presents difficulties when implementing it. In fact, determining the lower bound value of all coverings, as defined in Eq. 5, can be quite complex. For example, let's consider the uni-dimensional case, in which we want to compute the fd of a coastline (Koch Curve). According to Eq. 3 in the case of *D* = 1 the coastline length is measured by a ruler of length *r*. Accuracy of the measure increase with decreasing *r*. For *r* → 0 the coastline will have infinite length. Similar arguments can be applied to *<sup>D</sup>* <sup>=</sup> 2; for *<sup>r</sup>* <sup>→</sup> 0 the measure of *<sup>H</sup><sup>D</sup> <sup>δ</sup>* → 0.

This discussion implies that our coastline (ex. Koch Curve) will have a fd value more than one-dimensional and less than two- dimensional. For this reason, the fd is considered as the transition point (the lower bound value in Eq. 5) between *H<sup>D</sup> <sup>δ</sup>* <sup>→</sup> 0 and *<sup>H</sup><sup>D</sup> <sup>δ</sup>* → ∞.

Several computational approaches have been developed to avoid the need of defining the lower bound at issue. Therefore many strategies accomplished the fd computation by retrieving it from the scaling of the object's bulk with its size. In fact, object's bulk and its size have a linear relationship in a logarithmic scale so that the slope of the best fitting line may provide an accurate estimation of this relationship. By using this log-log graph, called *Richardson's plot*, the requirement of knowing the infimum over all coverings is relaxed.

Several approaches have been developed to estimate fractal dimension of images. In particular, this section will introduce two fractal analysis strategies: the *Box Counting Method* and the *Hand and Dividers Method*.

These methods overcome the problem by choosing as covering a simple rectangle fixed grid in order to obtain an upper bound on *DH*.

Five algorithms for a practical fd calculation based on these methods will also be presented.

#### **3.1. Box counting method**

The most popular method using the best fitting procedure is the so-called *Box Counting Method*[5][6]. Given a fractal structure *A* embedded in a d-dimensional volume the box-counting method basically consists of partitioning the structure space with a d-dimensional fixed-grid of square boxes of equal size *r*.

The number *N*(*r*) of nonempty boxes of size *r* needed to cover the fractal structure depends on *r*:

$$N(r) \sim r^{-D} \tag{8}$$

**Figure 2.** The Box-counting method applied to the Koch Curve with box size r = 0.4 (a); r = 1 (b); r = 1.4

The most popular methods are all based on the *Hand and Divers Method* which was originally

The Richardson method employs the so-called *walking technique* consisting of "walking"

In a nutshell, it corresponds to the length of the single step multiplied by the number of steps

The perimeter length of the boundary depends on the step length used so that a large step provides a rough estimation of the perimeter whereas a smaller step can take into account

In practice, the perimeter length is obtained by constructing a generally irregular polygon which approximate the border. Let *δB* be the set of coordinates of object boundary and let be

*l*(*�*) = *� n*(*�*) (10)

Fractal Dimension Estimation Methods for Biomedical Images 165

*Pi* = *l*(*�i*) = *�<sup>i</sup> n*(*�i*) (11)

*D* = 1 − *m* (12)

The actual structure boundary is so approximated by a polygon whose length is equal to:

introduced by Richardson[10] and successively developed by Mandelbrot[11].

around the boundary of the structure with a given step length.

The process is then reiterated for different step lengths:

With *Pi* the perimeter calculated with steps of length *�i*. The object's boundary fd *D* is finally estimate from:

Consequently, if the step length *�* decreases the perimeter *P* increases.

where *m* is the slope of the Richardson's plot.

(c); r= 2 (d)

needed to complete the walk.

finer details of the contour.

The box counting algorithm hence counts the number *N*(*r*) for different values of *r* and plot the log of the number *N*(*r*) versus the log of the actual box size *r*. The value of the box-counting dimension *D* is estimated from the Richardson's plot best fitting curve slope.

$$-D = \lim\_{r \to 0} \frac{\log N(r)}{\log r} \tag{9}$$

Figure 2 shows the Box counting method for the Koch Curve.

Several algorithms[7][8][9] based on box counting method have been developed and widely used for fd estimation, as it can be applied to sets with or without self-similarity. However, in computing fd with this method, one either counts or does not count a box according to whether there are no points or some points in the box. No provision is made for weighting the box according to the number of points belonging to the fractal and inside the current box.

### **3.2. Hand and dividers method**

Useful features and information can be deducted from the contours of structures belonging to an image and there is a number of techniques that can be used when estimating the boundary fractal dimension.

**Figure 2.** The Box-counting method applied to the Koch Curve with box size r = 0.4 (a); r = 1 (b); r = 1.4 (c); r= 2 (d)

The most popular methods are all based on the *Hand and Divers Method* which was originally introduced by Richardson[10] and successively developed by Mandelbrot[11].

The Richardson method employs the so-called *walking technique* consisting of "walking" around the boundary of the structure with a given step length.

The actual structure boundary is so approximated by a polygon whose length is equal to:

$$l(\epsilon) = \epsilon \ n(\epsilon) \tag{10}$$

In a nutshell, it corresponds to the length of the single step multiplied by the number of steps needed to complete the walk.

The process is then reiterated for different step lengths:

$$P\_{\mathbf{i}} = l(\boldsymbol{\varepsilon}\_{\mathbf{i}}) = \boldsymbol{\varepsilon}\_{\mathbf{i}} \ \boldsymbol{n}(\boldsymbol{\varepsilon}\_{\mathbf{i}}) \tag{11}$$

With *Pi* the perimeter calculated with steps of length *�i*.

The object's boundary fd *D* is finally estimate from:

$$D = 1 - m\tag{12}$$

where *m* is the slope of the Richardson's plot.

4 Will-be-set-by-IN-TECH

This discussion implies that our coastline (ex. Koch Curve) will have a fd value more than one-dimensional and less than two- dimensional. For this reason, the fd is considered as the

Several computational approaches have been developed to avoid the need of defining the lower bound at issue. Therefore many strategies accomplished the fd computation by retrieving it from the scaling of the object's bulk with its size. In fact, object's bulk and its size have a linear relationship in a logarithmic scale so that the slope of the best fitting line may provide an accurate estimation of this relationship. By using this log-log graph, called *Richardson's plot*, the requirement of knowing the infimum over all coverings is relaxed.

Several approaches have been developed to estimate fractal dimension of images. In particular, this section will introduce two fractal analysis strategies: the *Box Counting Method*

These methods overcome the problem by choosing as covering a simple rectangle fixed grid

Five algorithms for a practical fd calculation based on these methods will also be presented.

The most popular method using the best fitting procedure is the so-called *Box Counting Method*[5][6]. Given a fractal structure *A* embedded in a d-dimensional volume the box-counting method basically consists of partitioning the structure space with a

The number *N*(*r*) of nonempty boxes of size *r* needed to cover the fractal structure depends

The box counting algorithm hence counts the number *N*(*r*) for different values of *r* and plot the log of the number *N*(*r*) versus the log of the actual box size *r*. The value of the box-counting dimension *D* is estimated from the Richardson's plot best fitting curve slope.

Several algorithms[7][8][9] based on box counting method have been developed and widely used for fd estimation, as it can be applied to sets with or without self-similarity. However, in computing fd with this method, one either counts or does not count a box according to whether there are no points or some points in the box. No provision is made for weighting the box according to the number of points belonging to the fractal and inside the current box.

Useful features and information can be deducted from the contours of structures belonging to an image and there is a number of techniques that can be used when estimating the boundary

log *N*(*r*)

− *D* = lim *r*→0 *<sup>δ</sup>* <sup>→</sup> 0 and *<sup>H</sup><sup>D</sup>*

*<sup>N</sup>*(*r*) <sup>∼</sup> *<sup>r</sup>*−*<sup>D</sup>* (8)

log *<sup>r</sup>* (9)

*<sup>δ</sup>* → ∞.

transition point (the lower bound value in Eq. 5) between *H<sup>D</sup>*

and the *Hand and Dividers Method*.

**3.1. Box counting method**

**3.2. Hand and dividers method**

fractal dimension.

on *r*:

in order to obtain an upper bound on *DH*.

d-dimensional fixed-grid of square boxes of equal size *r*.

Figure 2 shows the Box counting method for the Koch Curve.

The perimeter length of the boundary depends on the step length used so that a large step provides a rough estimation of the perimeter whereas a smaller step can take into account finer details of the contour.

Consequently, if the step length *�* decreases the perimeter *P* increases.

In practice, the perimeter length is obtained by constructing a generally irregular polygon which approximate the border. Let *δB* be the set of coordinates of object boundary and let be  a fixed step length. Given a starting point, an arbitrary contour point (*xs*, *ys*), the next point on the boundary (*xs*<sup>2</sup> , *ys*<sup>2</sup> ) in a fixed direction (e.g. clockwise) is the point that has a distance

$$d\_l = \sqrt{(\mathbf{x}\_s - \mathbf{x}\_{s\_2})^2 + (y\_s - y\_{s\_2})^2} \tag{13}$$

**4.1. HYBRID algorithm**

4 shows the flow chart of the method.

the initial starting point *S* is reached.

The HYBRID algorithm is a computer implementation of Hand and Dividers method developed by Clark[12]. Let *δB* be the boundary of the object whose fd we wish to compute. The main part of the method focuses on the perimeter estimation and the corresponding Richardson's plot is then attained by reiterating this hard core part at different step size. Figure

Fractal Dimension Estimation Methods for Biomedical Images 167

**Figure 4.** Perimeter estimation by HYBRID method flowchart: an arbitrary *starting point S* (*xS*, *yS*) on the boundary line is chosen and copied in a new variable, which is called *current point C* (*xC*(*i*), *yC*(*i*)). The index i runs through the total number of coordinate points and is iteratively increased defining the *running point R* with coordinates (*xR*, *yR*). The distance d between *S* and *R* is calculated and a check on *d* when smaller than a fixed step length *�* is done. The process is repeated until a boundary point whose distance from (*xS*, *yS*) is larger than the step *�* is reached. The next pivot point on the boundary line is determined by choosing between the two points the closest to the step length. The distance is then stored and this point becomes the new starting point in order to calculate the next pivot point and so on, until

Given an arbitrary *starting point S* and its coordinates (*xS*, *yS*) on the boundary, the algorithm searches for the next pivot point. In particular the starting point is copied into a *current point*,

as close as possible to .

The reached point then becomes the new starting point and is used to locate the next point on the boundary that satisfies the previous condition. This process is repeated until the initial starting point is reached.

The sum of all distances *dj* corresponds to the irregular polygon perimeter (Figure3).

A number of different perimeters for each polygon at each fixed step length are used to build the Richardson's plot and the slope of its best linear fit is exploited to estimate the fd.

**Figure 3.** Walking technique applied to a coastline with different step lengths.

#### **4. Algorithms**

All Hand and dividers techniques rely on the same identical principle that attempt to approximate the border perimeter with a different polygons. However, since the point coordinates belonging to border set are discrete, all the implemented methods differ in the choice of which point in the set has a distance that better approximate the step length.

The following two methods are the implementations of two different choices about how to overcome this particular issue.

## **4.1. HYBRID algorithm**

6 Will-be-set-by-IN-TECH

 a fixed step length. Given a starting point, an arbitrary contour point (*xs*, *ys*), the next point on the boundary (*xs*<sup>2</sup> , *ys*<sup>2</sup> ) in a fixed direction (e.g. clockwise) is the point that has a distance

The reached point then becomes the new starting point and is used to locate the next point on the boundary that satisfies the previous condition. This process is repeated until the initial

A number of different perimeters for each polygon at each fixed step length are used to build

The sum of all distances *dj* corresponds to the irregular polygon perimeter (Figure3).

the Richardson's plot and the slope of its best linear fit is exploited to estimate the fd.

(*xs* − *xs*<sup>2</sup> )<sup>2</sup> + (*ys* − *ys*<sup>2</sup> )<sup>2</sup> (13)

*di* = 

**Figure 3.** Walking technique applied to a coastline with different step lengths.

All Hand and dividers techniques rely on the same identical principle that attempt to approximate the border perimeter with a different polygons. However, since the point coordinates belonging to border set are discrete, all the implemented methods differ in the choice of which point in the set has a distance that better approximate the step length.

The following two methods are the implementations of two different choices about how to

as close as possible to .

starting point is reached.

**4. Algorithms**

overcome this particular issue.

The HYBRID algorithm is a computer implementation of Hand and Dividers method developed by Clark[12]. Let *δB* be the boundary of the object whose fd we wish to compute. The main part of the method focuses on the perimeter estimation and the corresponding Richardson's plot is then attained by reiterating this hard core part at different step size. Figure 4 shows the flow chart of the method.

**Figure 4.** Perimeter estimation by HYBRID method flowchart: an arbitrary *starting point S* (*xS*, *yS*) on the boundary line is chosen and copied in a new variable, which is called *current point C* (*xC*(*i*), *yC*(*i*)). The index i runs through the total number of coordinate points and is iteratively increased defining the *running point R* with coordinates (*xR*, *yR*). The distance d between *S* and *R* is calculated and a check on *d* when smaller than a fixed step length *�* is done. The process is repeated until a boundary point whose distance from (*xS*, *yS*) is larger than the step *�* is reached. The next pivot point on the boundary line is determined by choosing between the two points the closest to the step length. The distance is then stored and this point becomes the new starting point in order to calculate the next pivot point and so on, until the initial starting point *S* is reached.

Given an arbitrary *starting point S* and its coordinates (*xS*, *yS*) on the boundary, the algorithm searches for the next pivot point. In particular the starting point is copied into a *current point*, *C* (*xC*, *yC*), which identifies all points having a mutual distance of about . The actual point running through the entire border is indicated as *running point R* (*xR*, *yR*).

point then becomes the new current point and is used to calculate the next boundary point

Fractal Dimension Estimation Methods for Biomedical Images 169

The process is stopped when we come back to the initial starting point (*xS*, *yS*) in order to

**Figure 5.** Perimeter estimation by EXACT method flowchart: an arbitrary *starting point S* (*xs*, *ys*) on the boundary line is chosen and copied in a new variable, which is called *current point C* (*xc* (*i*), *yc* (*i*)). The index i runs through the total number of coordinate points and is iteratively increased defining the *running point R* with coordinates (*xR*, *yR*). The distance d between *S* and *R* is calculated and a check on *d* when smaller than a fixed step length is done. The process is repeated until a boundary point whose distance from (*xs*, *ys*) is larger than the step is reached. The exact position of the next pivot point (*x*, *y*)

The point (*x*, *y*) becomes the new starting point in order to calculate the next pivot point and so on, until

The perimeter length of the polygon is found by adding the final incomplete step length to

The results, i.e., perimeter lengths versus step lengths, are plotted on a log-log Richardson's Plot. From the slope of the fitting line on the Richardson's plot we obtain the fd of the

on the boundary line is determined by interpolating the two consecutivepoints (*xR*, *yR*) and

the sum of the other step lengths needed to entirely cover the boundary.

The procedure is then repeated for different step lengths[15].

and so on.

(*xR*+1, *yR*<sup>+</sup>1).

the initial starting point *S* is reached.

examined boundary[1, 4, 12, 16–21]

obtain a polygon as is shown in Figure 8.

Therefore the program searches for a specific running point having a distance from *C* as near as possible to the step . In particular, in the HYBRID method the real step may be chosen to be longer or shorter than the fixed step depending on the minimum deviation from it. Similarly once the running point hits a contour point having a distance from the actual current one bigger than the size step, the choice is made between that point and the preceding one.

Afterward, the computed distance between these two points *R* and *C* is stored and the running point becomes the new current point.

The procedure continues until the initial starting point is reached. Obviously it is likely that after a complete walking the starting point *S* may be reached before having hit the following current point *C*. In other words, there may not be a multiple of step size so that the final incomplete step length *r* is added to the others stored distances, whose sum represent the boundary's perimeter. Since the fixed step length is adapted every time during the perimeter computation, its averaged value is then computed and used in the Richardson's plot.

## **4.2. EXACT algorithm**

The EXACT algorithm was proposed for the first time by Clark in 1986[12]. As it will be shown, this method requires a longer computational time by providing a simpler solution to the choice of the best current points.

Similarly to the previous method the entire perimeter estimation is displayed in the flow chart of Figure 5.

The procedure is very similar to the one used for the previous method. As before (see Figure 5), the end of the step may not coincide with the digitized coordinates of the boundary.

The way the EXACT method attempts to overcome this problem relies on the assumption of piecewise linearity, meaning that all the points on the contour can be joined by a series of straight line[13, 14] (see Figure 6 (a)).

The location of the next current point *C* on the boundary from the one previously determined is schematically illustrated in Figure 6 (b).

The procedure starts from an arbitrary starting point (*xS*, *yS*) and the algorithm searches for the next pivot point. In particular the starting point is copied into a *current point*, *C* (*xC*, *yC*), which identifies all points having a mutual distance of about . The actual point running through the entire border is indicated as *running point R* (*xR*, *yR*).

The distance from the current point to each point on the contour line is then calculated until the step length falls between two consecutive boundary points, (*xR*, *yR*) and (*xR*+1, *yR*+1) for which:

$$\sqrt{(\mathbf{x}\_{R} - \mathbf{x}\_{\mathbb{C}})^2 + (y\_{R} - y\_{\mathbb{C}})^2} < \epsilon < \sqrt{(\mathbf{x}\_{R+1} - \mathbf{x}\_{\mathbb{C}})^2 + (y\_{R+1} - y\_{\mathbb{C}})^2} \tag{14}$$

The exact position of the point *N* with coordinates (*x*, *y*) is deduced by a process of geometric interpolation between the two consecutive running points (*xR*, *yR*) and (*xR*+1, *yR*+1). This point then becomes the new current point and is used to calculate the next boundary point and so on.

8 Will-be-set-by-IN-TECH

*C* (*xC*, *yC*), which identifies all points having a mutual distance of about . The actual point

Therefore the program searches for a specific running point having a distance from *C* as near as possible to the step . In particular, in the HYBRID method the real step may be chosen to be longer or shorter than the fixed step depending on the minimum deviation from it. Similarly once the running point hits a contour point having a distance from the actual current one bigger than the size step, the choice is made between that point and the preceding one.

Afterward, the computed distance between these two points *R* and *C* is stored and the running

The procedure continues until the initial starting point is reached. Obviously it is likely that after a complete walking the starting point *S* may be reached before having hit the following current point *C*. In other words, there may not be a multiple of step size so that the final incomplete step length *r* is added to the others stored distances, whose sum represent the boundary's perimeter. Since the fixed step length is adapted every time during the perimeter

The EXACT algorithm was proposed for the first time by Clark in 1986[12]. As it will be shown, this method requires a longer computational time by providing a simpler solution to

Similarly to the previous method the entire perimeter estimation is displayed in the flow chart

The procedure is very similar to the one used for the previous method. As before (see Figure 5), the end of the step may not coincide with the digitized coordinates of the boundary.

The way the EXACT method attempts to overcome this problem relies on the assumption of piecewise linearity, meaning that all the points on the contour can be joined by a series of

The location of the next current point *C* on the boundary from the one previously determined

The procedure starts from an arbitrary starting point (*xS*, *yS*) and the algorithm searches for the next pivot point. In particular the starting point is copied into a *current point*, *C* (*xC*, *yC*), which identifies all points having a mutual distance of about . The actual point running

The distance from the current point to each point on the contour line is then calculated until the step length falls between two consecutive boundary points, (*xR*, *yR*) and (*xR*+1, *yR*+1)

The exact position of the point *N* with coordinates (*x*, *y*) is deduced by a process of geometric interpolation between the two consecutive running points (*xR*, *yR*) and (*xR*+1, *yR*+1). This

(*xR*+<sup>1</sup> − *xC*)<sup>2</sup> + (*yR*+<sup>1</sup> − *yC*)<sup>2</sup> (14)

computation, its averaged value is then computed and used in the Richardson's plot.

running through the entire border is indicated as *running point R* (*xR*, *yR*).

point becomes the new current point.

**4.2. EXACT algorithm**

of Figure 5.

for which:

the choice of the best current points.

straight line[13, 14] (see Figure 6 (a)).

is schematically illustrated in Figure 6 (b).

through the entire border is indicated as *running point R* (*xR*, *yR*).

(*xR* − *xC*)<sup>2</sup> + (*yR* − *yC*)<sup>2</sup> < <

The process is stopped when we come back to the initial starting point (*xS*, *yS*) in order to obtain a polygon as is shown in Figure 8.

**Figure 5.** Perimeter estimation by EXACT method flowchart: an arbitrary *starting point S* (*xs*, *ys*) on the boundary line is chosen and copied in a new variable, which is called *current point C* (*xc* (*i*), *yc* (*i*)). The index i runs through the total number of coordinate points and is iteratively increased defining the *running point R* with coordinates (*xR*, *yR*). The distance d between *S* and *R* is calculated and a check on *d* when smaller than a fixed step length is done. The process is repeated until a boundary point whose distance from (*xs*, *ys*) is larger than the step is reached. The exact position of the next pivot point (*x*, *y*) on the boundary line is determined by interpolating the two consecutivepoints (*xR*, *yR*) and (*xR*+1, *yR*<sup>+</sup>1).

The point (*x*, *y*) becomes the new starting point in order to calculate the next pivot point and so on, until the initial starting point *S* is reached.

The perimeter length of the polygon is found by adding the final incomplete step length to the sum of the other step lengths needed to entirely cover the boundary.

The procedure is then repeated for different step lengths[15].

The results, i.e., perimeter lengths versus step lengths, are plotted on a log-log Richardson's Plot. From the slope of the fitting line on the Richardson's plot we obtain the fd of the examined boundary[1, 4, 12, 16–21]

10 Will-be-set-by-IN-TECH 170 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 3 Fractal Dimension Estimation Methods for Biomedical Images <sup>11</sup>

**Figure 8.** The Box-counting algorithm flowchart: given an image *I*, its size, *S*0, is set as the maximum size from which the computer program starts to calculate the others decreasing box-sizes according to *S* = *S*0/2*<sup>N</sup>*. The *S* value has minimum value which is equal to the pixel size. The number *p* indicates the total number of box size *S*. The next step is a check on whether at least one pixel is in the box: if the box

Fractal Dimension Estimation Methods for Biomedical Images 171

maximum index *p*(*s*) is reached. Then the number *Nboxes*(*S*) for a given size *S* is stored and the process restarts with a different box size. When the minimum box size is reached the program stops and gives the output variables of *Nboxes* and the size value. Using the Eq. 8 the fractal dimension *D* can be

Figure 8 shows the flow chart for box-counting fd estimation and for different box sizes. Moreover, since the procedure of size scaling (*S* = *S*0/2*<sup>N</sup>* with *N* number of iterations) may be not always applicable to any image matrix size, image padding with background pixels is

Therefore the final image *I* has a dimension that is a power of 2. This can be easily

The box counting method is an extremely powerful tool for fd computation; in fact, it is easy

However, a major limitation lies on the fact that the counting process of nonempty boxes implies its use only for binary images rather than gray scale ones. An extension of the

is non-empty, the check is stopped when one pixel is found. The procedure continues until the

estimate, from the least square linear fit.

implemented by using *padarray* matlab function.

to implement as well as flexible and robust.

**4.4. Differential Box-counting algorithm (DBC)**

performed.

**Figure 6.** a) The piece-wise linear assumption (a) (b) and the EXACT algorithm (c); b)Geometric EXACT interpolation scheme, with *S* starting point given by (*xC*, *yC*) coordinates, *R* and *R* two consecutive boundary running points respectly given by (*xR*, *yR*) and (*xR*+1, *yR*<sup>+</sup>1) coordinates, *N* new current point obtained by the interpolation between *R* and *R* and a fixed step length.

**Figure 7.** MRI image of an aneurismatic bone cyst (a), (b). Walking technique applied to an aneurysmatic bone cyst boundary (c).

#### **4.3. Box-counting algorithm**

The Box-counting algorithm implementation of box-counting method relies on the basic idea of covering a given digital binary image with a set of measuring boxes of sizes *S* and then to count the number of boxes which actually contain the image.

10 Will-be-set-by-IN-TECH

(a) (b)

**Figure 6.** a) The piece-wise linear assumption (a) (b) and the EXACT algorithm (c); b)Geometric EXACT interpolation scheme, with *S* starting point given by (*xC*, *yC*) coordinates, *R* and *R* two consecutive boundary running points respectly given by (*xR*, *yR*) and (*xR*+1, *yR*<sup>+</sup>1) coordinates, *N* new current point

**Figure 7.** MRI image of an aneurismatic bone cyst (a), (b). Walking technique applied to an

count the number of boxes which actually contain the image.

The Box-counting algorithm implementation of box-counting method relies on the basic idea of covering a given digital binary image with a set of measuring boxes of sizes *S* and then to

aneurysmatic bone cyst boundary (c).

**4.3. Box-counting algorithm**

obtained by the interpolation between *R* and *R* and a fixed step length.

**Figure 8.** The Box-counting algorithm flowchart: given an image *I*, its size, *S*0, is set as the maximum size from which the computer program starts to calculate the others decreasing box-sizes according to *S* = *S*0/2*<sup>N</sup>*. The *S* value has minimum value which is equal to the pixel size. The number *p* indicates the total number of box size *S*. The next step is a check on whether at least one pixel is in the box: if the box is non-empty, the check is stopped when one pixel is found. The procedure continues until the maximum index *p*(*s*) is reached. Then the number *Nboxes*(*S*) for a given size *S* is stored and the process restarts with a different box size. When the minimum box size is reached the program stops and gives the output variables of *Nboxes* and the size value. Using the Eq. 8 the fractal dimension *D* can be estimate, from the least square linear fit.

Figure 8 shows the flow chart for box-counting fd estimation and for different box sizes. Moreover, since the procedure of size scaling (*S* = *S*0/2*<sup>N</sup>* with *N* number of iterations) may be not always applicable to any image matrix size, image padding with background pixels is performed.

Therefore the final image *I* has a dimension that is a power of 2. This can be easily implemented by using *padarray* matlab function.

### **4.4. Differential Box-counting algorithm (DBC)**

The box counting method is an extremely powerful tool for fd computation; in fact, it is easy to implement as well as flexible and robust.

However, a major limitation lies on the fact that the counting process of nonempty boxes implies its use only for binary images rather than gray scale ones. An extension of the standard approach to gray scale images is called the *Differential Box Counting (DBC)* and has been proposed in 1994 by N. Sarkar and Chaudhuri[8].

In the DBC method, a gray level image *I* is considered as a 3-D spatial surface with (*x*, *y*) denoting the pixels spatial coordinates and the third axis *z* the pixels gray level.

As for the standard box counting, the *M* × *M* image matrix is partitioned into non-overlapping *s* × *s*-sized boxes, where *s* is an integer falling in the interval [*M*/2, 1].

Then, the scale of each block is *r* = *s*. On each block there is a column of boxes of size *s* × *s* × *s*� , with *s*� denoting the height of a single box. Named *G* the total number of gray levels in *I*, hence *s*� is defined by the relationship *G*/*s*� = *M*/*s*[7].

Let numbers 1, 2, 3... be assigned to the boxes so to group the gray levels. Let the minimum and the maximum gray level of the image in the (*i*, *j*)*th* grid fall in box number *k* and *l*, respectively.

The number of boxes covering this block is calculated as:

$$n\_r(i,j) = l - k + 1\tag{15}$$

The DBC procedure has some weak points in the method used to select an appropriate box height[7], since the values of *s* is limited to the image size and *s*� is limited by the number of

Secondly, the box number calculation may lead to overestimate the number of boxes needed to cover the surface. Let *A* and *B* be the pixels associated with the minimum and the maximum

**Figure 10.** Example of DBC method application with boxes of *s* × *s* × *s*, when *s* = 3. The two pixels *A* and *B*, denoting the maximum and the minimum gray levels of the block, are assigned in two differents

According to DBC procedure, the two pixels are assigned in boxes 2 and boxes 3. The distance

Hence, when calculating Eq. 15, the block can be covered by a single box but its pixels with

To solve the aforementioned problems some modifications was proposed by J. Li, Q. Du and C. Sun[7]. Given a digital image *I* of size *M* × *M*, a new scale *r* is defined instead of *r*, i.e.

In particular, let *μ* and *σ* be the mean and the standard deviation of *I* respectively. Hence, if the greater part of image pixels fall into the interval of gray level within [*μ* − *aσ*, *μ* + *aσ*],

*<sup>r</sup>*� <sup>=</sup> *<sup>r</sup>*

If *dz* is the height of the boxes in the direction of *z*, the number of the column of boxes on a single image block correspond to the integer part of (*dz*/*r*� + 1) instead of (*dz*/*r* + 1) as in the original DBC method. Thus, since *r*� < *r*, the residual part of *dz*/*r*� is smaller than that of *dz*/*r*. As a result, the errors introduced using *r*� are smaller than in the original DBC method. A box with smaller height is chosen when a higher intensity variation is present on the image

<sup>1</sup> <sup>+</sup> <sup>2</sup>*a<sup>σ</sup>* (17)

Fractal Dimension Estimation Methods for Biomedical Images 173

blocks of size *s* × *s* in which the image is divided.

gray level of the block respectively, as is illustrated in Figure 10.

boxes, having distance in eight direction smaller than the box size *s* = 3.

between *A* and *B* is smaller than 3, which is the size of the box.

minimum and maximum gray levels fall into two different boxes.

where *a* is a positive integer, the height of the boxes is given by:

surface. So the improved method uses, in general, finer scales to count[7].

*r*� = *r*/*c* where *c*1 is a positive real number.

In Figure 9 for example *s* = *s*� = 3, hence *nr* = 3 − 1 + 1. Extending to the contribution from

all blocks:

$$N\_r = \sum\_{i,j} n\_r(i,j) \tag{16}$$

The Eq. 16 is computed for different box size *s* (so for different *r*) and the values of *Nr* are plotted versus the values of *r* in a log-log plot.

A matlab implementation of DBC can make use of functions such as *blockproc* or *col filt* in order to make the box partitioning and apply the Eq. 15.

The DBC procedure has some weak points in the method used to select an appropriate box height[7], since the values of *s* is limited to the image size and *s*� is limited by the number of blocks of size *s* × *s* in which the image is divided.

12 Will-be-set-by-IN-TECH

standard approach to gray scale images is called the *Differential Box Counting (DBC)* and has

In the DBC method, a gray level image *I* is considered as a 3-D spatial surface with (*x*, *y*)

As for the standard box counting, the *M* × *M* image matrix is partitioned into non-overlapping

Then, the scale of each block is *r* = *s*. On each block there is a column of boxes of size *s* × *s* × *s*�

with *s*� denoting the height of a single box. Named *G* the total number of gray levels in *I*, hence

Let numbers 1, 2, 3... be assigned to the boxes so to group the gray levels. Let the minimum and the maximum gray level of the image in the (*i*, *j*)*th* grid fall in box number *k* and *l*,

In Figure 9 for example *s* = *s*� = 3, hence *nr* = 3 − 1 + 1. Extending to the contribution from

**Figure 9.** Example of DBC method application for determining the number of boxes of size *s* × *s* × *s*,

*Nr* = ∑ *i*,*j*

The Eq. 16 is computed for different box size *s* (so for different *r*) and the values of *Nr* are

A matlab implementation of DBC can make use of functions such as *blockproc* or *col filt* in

*nr*(*i*, *j*) = *l* − *k* + 1 (15)

*nr*(*i*, *j*) (16)

,

denoting the pixels spatial coordinates and the third axis *z* the pixels gray level.

*s* × *s*-sized boxes, where *s* is an integer falling in the interval [*M*/2, 1].

been proposed in 1994 by N. Sarkar and Chaudhuri[8].

*s*� is defined by the relationship *G*/*s*� = *M*/*s*[7].

plotted versus the values of *r* in a log-log plot.

order to make the box partitioning and apply the Eq. 15.

The number of boxes covering this block is calculated as:

respectively.

when *s* = 3. all blocks:

Secondly, the box number calculation may lead to overestimate the number of boxes needed to cover the surface. Let *A* and *B* be the pixels associated with the minimum and the maximum gray level of the block respectively, as is illustrated in Figure 10.

**Figure 10.** Example of DBC method application with boxes of *s* × *s* × *s*, when *s* = 3. The two pixels *A* and *B*, denoting the maximum and the minimum gray levels of the block, are assigned in two differents boxes, having distance in eight direction smaller than the box size *s* = 3.

According to DBC procedure, the two pixels are assigned in boxes 2 and boxes 3. The distance between *A* and *B* is smaller than 3, which is the size of the box.

Hence, when calculating Eq. 15, the block can be covered by a single box but its pixels with minimum and maximum gray levels fall into two different boxes.

To solve the aforementioned problems some modifications was proposed by J. Li, Q. Du and C. Sun[7]. Given a digital image *I* of size *M* × *M*, a new scale *r* is defined instead of *r*, i.e. *r*� = *r*/*c* where *c*1 is a positive real number.

In particular, let *μ* and *σ* be the mean and the standard deviation of *I* respectively. Hence, if the greater part of image pixels fall into the interval of gray level within [*μ* − *aσ*, *μ* + *aσ*], where *a* is a positive integer, the height of the boxes is given by:

$$r' = \frac{r}{1 + 2a\sigma} \tag{17}$$

If *dz* is the height of the boxes in the direction of *z*, the number of the column of boxes on a single image block correspond to the integer part of (*dz*/*r*� + 1) instead of (*dz*/*r* + 1) as in the original DBC method. Thus, since *r*� < *r*, the residual part of *dz*/*r*� is smaller than that of *dz*/*r*.

As a result, the errors introduced using *r*� are smaller than in the original DBC method. A box with smaller height is chosen when a higher intensity variation is present on the image surface. So the improved method uses, in general, finer scales to count[7].

#### 14 Will-be-set-by-IN-TECH 174 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 3 Fractal Dimension Estimation Methods for Biomedical Images <sup>15</sup>

Moreover, the use of *dz* instead of *z* to count the number of boxes leads to the following modification of Eq. 15:

$$m\_r = \begin{cases} cei l(\frac{l-k}{r'}), \\ 1, & l=k \end{cases} \tag{18}$$

On the table 2 the computation results for the box counting method are also shown. The type of the displayed values are similar to the previous ones with the exception of Box counting uncertainty. In fact, the way an image can be partitioned into several boxes may affect the

Fractal Dimension Estimation Methods for Biomedical Images 175

To investigate the variability of the fd for different box partitioning layouts, random box subdivisions have been applied. Therefore, the results on the table 2 show the standard deviation of the different computed fds and the mean values for each fractal at issue. In

*Fractal f dtheo f dexp Time* (sec) *BC error Image size* Apollonian Gasket 1.3057 1.408 1.5 0.001 2000 × 2000 Sierpinski 1.5849 1.587 0.3 0.005 1000 × 1000 Dragon 2.0000 1.747 7.2 0.006 3670 × 3978 Hexaflake 1.7719 1.640 1.6 0.011 1050 × 1050

*Fractal f dtheo f dexp Time* (sec) *BC error Vector size*

Twin Dragon Hybrid 1.5236 1.466 8.6 0.006 117005 Twin Dragon Exact 1.5236 1.465 11.5 0.006 117005 Dragon Hybrid 1.5236 1.474 11.1 0.005 115665 Dragon Exact 1.5236 1.462 12.8 0.004 115665 Koch Hybrid 1.2619 1.276 31.2 0.004 786433 Koch Exact 1.2619 1.260 154.9 0.003 786433 Gosper Hybrid 1.1292 1.133 3.8 0.001 23280 Gosper Exact 1.1292 1.128 4.7 0.001 23280

In general, the EXACT and the HYBRID methods appeared to be more precise than the box counting method but on the other hand they have a less wide range of applicability. However, this is also the reason of the fortune of the box counting methods compared to the others. Also, HYBRID technique is computationally less expensive than EXACT especially when the number of border points is quite large. The use of a variable step length which can be shorter or longer than the fixed step size leads to a larger variability and so to a Richardoson's plot having a less accurate fitting. That has effects on the uncertainties of the parameter to estimate. Because of that, a more careful choice of the step size range is needed in the case of HYBRID

Importantly, it is quite clear that the choice of the starting point may also affect the perimeter value as the following currents points will depend upon this. A test on 80 random starting points for the Gosper Island fractal revealed that the fd computation performed with the

As for walking method, in box counting the process of scaling from the maximum box size is limited by the pixel size so in principle a gross resolution might be the reason of a bad estimate of fd. It is noteworthy that the tests performed do not show any correlation between resolution and fd accuracy; that may be also caused by the fact that some fractals such as

HYBRID method appeared to be more stable than the one with EXACT.

dragon does not reproduce the real fractal at small scales.

general, that variability is more pronounced in images having rougher resolution.

final computation of the number of nonempty boxes.

**Table 1.** Tabular of results for box counting method application.

**Table 2.** Tabular of results for walking-based methods application.

method.

with ceil( . ) denoting the function rounds the elements of the quantity into ( . ) to the nearest integers greater or equal to it.

Eq. 18 relies a new way to count the number of boxes that cover the (*i*, *j*)*th* block surface in which the boxes are assigned to the minimum gray level to the block rather than gray level 0[7].

As an example, suppose that the (*i*, *j*)*th* block is covered by a column boxes with the size 3*x*3*x*3. If the pixels *A* and *B* represent the maximum and the minimum gray levels of the block, the two pixels will be assigned as in Figure 10.

According to Eq. 18 the number of counted boxes is *nr* = 1, which is exactly the number of boxes covering the block.

As in standard box counting method, after having determined the number nr(i,j) for each block, the total number of boxes *Nr* covering the full image surface is computed for different scales *r*. Plotting the linear fit of log *Nr* versus the log *r* (Richardson's plot) the fd is finally estimated.

## **5. Applications and discussion**

Each described method has been implemented in Matlab 2010a and applied to either well-known fractals or biomedical images.

The results on the hand and dividers methods are shown in the table 1. The computed values are also compared to the theoretical fd values. The computational time for a 2.50 GHz 5i CPU is also shown.

The value ranges for the step size are not displayed but they were automatically chosen based upon the computation of the structure's maximum caliber diameter which is defined as the major axis of an ellipse in which the structure can be embedded. The range was then running from the 40% of the maximum caliber diameter to the minimum step defined as the maximum distance between any two contiguous border points.

In practice, both EXACT and HYBRID methods computed the different step sizes by scaling each time the maximum step by *a<sup>k</sup>* with *k* the number of the iteration. The chosen value of *a* = 1.2 is a compromise between a sufficient number of fitting points and the need to avoid too small variations of the step size so to duplicate perimeter estimation. The latter usually occurs in HYBRID method for it hits the same current points if the step does not vary enough in two consecutive iterations.

The parameter's estimation uncertainty is also shown in the table 1; that is calculated from the fitting accuracy based upon standard linear regression.

The number of data points used in the Richardson's plot was about 60 and two examples of that computation using EXACT and HYBRID are shown in Figure 12.

On the table 2 the computation results for the box counting method are also shown. The type of the displayed values are similar to the previous ones with the exception of Box counting uncertainty. In fact, the way an image can be partitioned into several boxes may affect the final computation of the number of nonempty boxes.

To investigate the variability of the fd for different box partitioning layouts, random box subdivisions have been applied. Therefore, the results on the table 2 show the standard deviation of the different computed fds and the mean values for each fractal at issue. In general, that variability is more pronounced in images having rougher resolution.



**Table 1.** Tabular of results for box counting method application.

14 Will-be-set-by-IN-TECH

Moreover, the use of *dz* instead of *z* to count the number of boxes leads to the following

with ceil( . ) denoting the function rounds the elements of the quantity into ( . ) to the nearest

Eq. 18 relies a new way to count the number of boxes that cover the (*i*, *j*)*th* block surface in which the boxes are assigned to the minimum gray level to the block rather than gray level

As an example, suppose that the (*i*, *j*)*th* block is covered by a column boxes with the size 3*x*3*x*3. If the pixels *A* and *B* represent the maximum and the minimum gray levels of the

According to Eq. 18 the number of counted boxes is *nr* = 1, which is exactly the number of

As in standard box counting method, after having determined the number nr(i,j) for each block, the total number of boxes *Nr* covering the full image surface is computed for different scales *r*. Plotting the linear fit of log *Nr* versus the log *r* (Richardson's plot) the fd is finally

Each described method has been implemented in Matlab 2010a and applied to either

The results on the hand and dividers methods are shown in the table 1. The computed values are also compared to the theoretical fd values. The computational time for a 2.50 GHz 5i CPU

The value ranges for the step size are not displayed but they were automatically chosen based upon the computation of the structure's maximum caliber diameter which is defined as the major axis of an ellipse in which the structure can be embedded. The range was then running from the 40% of the maximum caliber diameter to the minimum step defined as the maximum

In practice, both EXACT and HYBRID methods computed the different step sizes by scaling each time the maximum step by *a<sup>k</sup>* with *k* the number of the iteration. The chosen value of *a* = 1.2 is a compromise between a sufficient number of fitting points and the need to avoid too small variations of the step size so to duplicate perimeter estimation. The latter usually occurs in HYBRID method for it hits the same current points if the step does not vary enough

The parameter's estimation uncertainty is also shown in the table 1; that is calculated from the

The number of data points used in the Richardson's plot was about 60 and two examples of

1, *<sup>l</sup>* <sup>=</sup> *<sup>k</sup>* (18)

*ceil*(*l*−*<sup>k</sup> <sup>r</sup>*� ),

*nr* =

block, the two pixels will be assigned as in Figure 10.

modification of Eq. 15:

0[7].

estimated.

is also shown.

in two consecutive iterations.

integers greater or equal to it.

boxes covering the block.

**5. Applications and discussion**

well-known fractals or biomedical images.

distance between any two contiguous border points.

fitting accuracy based upon standard linear regression.

that computation using EXACT and HYBRID are shown in Figure 12.

**Table 2.** Tabular of results for walking-based methods application.

In general, the EXACT and the HYBRID methods appeared to be more precise than the box counting method but on the other hand they have a less wide range of applicability. However, this is also the reason of the fortune of the box counting methods compared to the others. Also, HYBRID technique is computationally less expensive than EXACT especially when the number of border points is quite large. The use of a variable step length which can be shorter or longer than the fixed step size leads to a larger variability and so to a Richardoson's plot having a less accurate fitting. That has effects on the uncertainties of the parameter to estimate. Because of that, a more careful choice of the step size range is needed in the case of HYBRID method.

Importantly, it is quite clear that the choice of the starting point may also affect the perimeter value as the following currents points will depend upon this. A test on 80 random starting points for the Gosper Island fractal revealed that the fd computation performed with the HYBRID method appeared to be more stable than the one with EXACT.

As for walking method, in box counting the process of scaling from the maximum box size is limited by the pixel size so in principle a gross resolution might be the reason of a bad estimate of fd. It is noteworthy that the tests performed do not show any correlation between resolution and fd accuracy; that may be also caused by the fact that some fractals such as dragon does not reproduce the real fractal at small scales.

(a) (b)

Fractal Dimension Estimation Methods for Biomedical Images 177

**Figure 13.** High resolution mammography image (a); fd recostruction image by standard Differential

In this chapter some of the most widely used and robust methods for fractal dimension estimation as well as their performances have been described. For few of them a detailed description of the algorithm has been also reported to make much easier for a beginner to start and implement his own Matlab code. Computational time is not excessively long to necessitate compiled functions such as C-mex files but that can be an advantage when using very high resolution images. The use of the described algorithms is obviously not restricted to the sole field of the image processing but it can be applied with some changes to any data

(c)

**6. Conclusions**

analysis.

Box Counting (DBC) (b); fd reconstruction image by modified DBC (c).

**Figure 11.** EXACT method apllied to the twin dragon fractal: Richardson's Plot.

**Figure 12.** HYBRID method apllied to the twin dragon fractal: Richardson's Plot.

An application of the DBC method on a x-ray image is also shown in Figure 13 where breast cancer mammography image has been processed. The method uses a sliding technique as implemented in *blockproc* or *col filt* matlab functions so to produce an image rather than a single fd value as previously described.

The second DBC method shows higher contrast in the area of the cancer and consequently lower fd values. Due to the enormous amount of linear fitting performed for an image size of 3450 × 3100 the computational time reached 15 minutes.

16 Will-be-set-by-IN-TECH

**Figure 11.** EXACT method apllied to the twin dragon fractal: Richardson's Plot.

**Figure 12.** HYBRID method apllied to the twin dragon fractal: Richardson's Plot.

single fd value as previously described.

3450 × 3100 the computational time reached 15 minutes.

An application of the DBC method on a x-ray image is also shown in Figure 13 where breast cancer mammography image has been processed. The method uses a sliding technique as implemented in *blockproc* or *col filt* matlab functions so to produce an image rather than a

The second DBC method shows higher contrast in the area of the cancer and consequently lower fd values. Due to the enormous amount of linear fitting performed for an image size of

(a) (b)

(c)

**Figure 13.** High resolution mammography image (a); fd recostruction image by standard Differential Box Counting (DBC) (b); fd reconstruction image by modified DBC (c).

### **6. Conclusions**

In this chapter some of the most widely used and robust methods for fractal dimension estimation as well as their performances have been described. For few of them a detailed description of the algorithm has been also reported to make much easier for a beginner to start and implement his own Matlab code. Computational time is not excessively long to necessitate compiled functions such as C-mex files but that can be an advantage when using very high resolution images. The use of the described algorithms is obviously not restricted to the sole field of the image processing but it can be applied with some changes to any data analysis.

## **Author details**

Antonio Napolitano, Sara Ungania and Vittorio Cannata

*Department of Occupational Health and Safety, Medical Physics, Bambino Gesù Children's Hospital, Rome, Italy*

**MATLAB Aided Option Replication**

In reality markets are incomplete in the sense that perfect replication of contingent claims using only the underlying asset and a riskless bond is impossible. In other words, that is perfect risk transfer is not possible since some payoffs cannot be replicated by trading in marketed securities. From the work of Ross in [21], it is evident that whenever the payoff of every call or put option can be replicated then the securities market is complete. In addition, an important implication of the aforementioned work of Ross, is the existence of options that cannot be replicated by the primitive securities when markets are incomplete. In [7], the authors came to the conclusion that Ross's result is, in fact, a negative result since it asserts that in an incomplete market one cannot expect to replicate the payoff of each option even if the underlying asset is traded. In the same paper, it is proved that if the number of securities is less than half the number of states of the world, then (generically) not a single option can be replicated by traded securities. In [10], the author extended the aforementioned result in [7], to accommodate cases where the condition on the number of primitive securities is not imposed. In particular, it is proved that if there exists no binary payoff vector in the asset span, then for each portfolio there exists a set of nontrivial exercise prices of full measure such that any option on the portfolio with an exercise price in this set is non-replicated. Furthermore, note that the results of Ross for two-date security markets with finitely many states holds for

**Chapter 8**

It is well accepted that the lattice theoretic ideas are the most important technical contributions of the large literature on infinite-dimensional modern mathematical finance (for example lattice theoretic ideas in general equilibrium theory). However, ordered vector spaces that are not lattice ordered arise naturally in models of portfolio trading. Moreover, if available securities have smooth payoffs, then the portfolio space is never a vector lattice. It should be pointed out that since call and put options are vector lattice operations in the space of contingent claims, their replication by available securities requires a vector lattice structure in the portfolio space. There is a large literature on vector lattice theory related to mathematical

On the other hand, there is an obvious need for properly structured high performance computational methods on vector lattices. Moreover, the main concern is to describe, in

> ©2012 Katsikis, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly

©2012 Katsikis, licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Additional information is available at the end of the chapter

security markets with more than two dates, see [8, 9].

economics; see for instance [1–7, 17–20, 22, 23].

cited.

Vasilios N. Katsikis

**1. Introduction**

http://dx.doi.org/10.5772/45723

## **7. References**


## **MATLAB Aided Option Replication**

Vasilios N. Katsikis

18 Will-be-set-by-IN-TECH

*Department of Occupational Health and Safety, Medical Physics, Bambino Gesù Children's Hospital,*

[1] J. Orford and W. Whalley, *The use of the fractal dimension to quantify the morphology of*

[2] J. Keller et al., *Texture description and segmentation through fractal geometry*, [Computer

[6] S. Deepa,T. Tessamma *Fractal Features based on Differential Box Counting Method for the Categoritazion of Digital Mammograms*, [International Journal of Computer Information

[7] J. Li, Q. Du, *An improved box-counting method for image fractal dimension estimation*, [

[8] N. Sarker, B. B. Chaudhuri, *An efficient differential box-counting approach to compute fractal dimension of image*, [IEEE Transaction on Systems, Man, and Cybernetics, 24, 115-120 ],

[9] A. P. Pentland, *Fractal-based description of natural scenes*, [IEEE Transaction on Pattern

[10] L. F. Richardson, *Fractal growth phenomena*, [Ann. Arbor, Mich. : The Society 6, 139] (1961). [11] B. Mandelbrot, *How long is the coast of Britain? Statistical self-similarity and fractional*

[12] N. Clark, *Three techiques for implementing digital fractal analysis of particle shape*, [Powder

[13] M. Allen, G. J. Brown, N. J. Miles, *Measurement of boundary fractal dimensions: review of*

[15] P. Podsiadlo, G. W. Stachowiak, *Evaluation of boundary fractal methods for the*

[16] J. D. Farmer, E. Ott. and J. A. Yorke, *The dimension of chaotic attractors*, [Physica 7D,

[19] H. Schwarz and E. Exner, *The implementation of the concept of fractal dimension on*

[21] H. von Koch, *Sur une courbe continue sans tangente, obtenue par une construction geometrique*

*semi-automatic image analyzer*, [Powder Technology, 27, 207-213], (1980).

*elementaire*, Archiv for Matemat., [Astron. och Fys. 1, 681-702] (1904).

[17] B. Kaye, *A Random Walk Through Fractal Dimension*, [VCH Verlagsgesellschaft] (1986) [18] L. Niemeyer, L. Pietronero and H. J. Wiesmann, *Fractals dimension of dielectric breakdown*,

*irregular-shaped particles*, [Sedimentology, 30, 655-668], (1983).

Vision, Graphics and Image Processing, 45, 150-166 ], (1989).

[3] F. Hausdorff, *Dimension und ausseres Mass*, [Math. Annalen 79, 157] (1919). [4] J. Theiler, *Estimating fractal dimension*, [7, 6/June 1990/J. Opt. Soc. Am. A] (1989). [5] B. Mandelbrot, *The Fractal Geometry of Nature*, W.H.Freeman and Company, New York,

System and Industrial Management Applications, 2, 011-019 ], ( 2010).

**Author details**

**7. References**

(1983 ).

(1994).

153] (1983).

*Rome, Italy*

Antonio Napolitano, Sara Ungania and Vittorio Cannata

Pattern Recognition, 42, 2460-2469 ], (2009).

*dimension*, [Science 155, 636-638] (1967).

Technology 46, 45] (1986).

[Phys. Rev. Lett. 52, 1033] (1984).

[20] J. Russ, *Fractals surfaces*, [Plenum Press] (1994).

Analysis and Machine Intelligence, 6, 661-674 ], (1984).

*current techinques*, University of Nottingham , UK (1994). [14] M. Allen, *Ph.D. Thesis*, University of Nottingham , UK (1994).

*characterization of wear particles*, [Wear 217(1) , 24-34] (1998).

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/45723

## **1. Introduction**

In reality markets are incomplete in the sense that perfect replication of contingent claims using only the underlying asset and a riskless bond is impossible. In other words, that is perfect risk transfer is not possible since some payoffs cannot be replicated by trading in marketed securities. From the work of Ross in [21], it is evident that whenever the payoff of every call or put option can be replicated then the securities market is complete. In addition, an important implication of the aforementioned work of Ross, is the existence of options that cannot be replicated by the primitive securities when markets are incomplete. In [7], the authors came to the conclusion that Ross's result is, in fact, a negative result since it asserts that in an incomplete market one cannot expect to replicate the payoff of each option even if the underlying asset is traded. In the same paper, it is proved that if the number of securities is less than half the number of states of the world, then (generically) not a single option can be replicated by traded securities. In [10], the author extended the aforementioned result in [7], to accommodate cases where the condition on the number of primitive securities is not imposed. In particular, it is proved that if there exists no binary payoff vector in the asset span, then for each portfolio there exists a set of nontrivial exercise prices of full measure such that any option on the portfolio with an exercise price in this set is non-replicated. Furthermore, note that the results of Ross for two-date security markets with finitely many states holds for security markets with more than two dates, see [8, 9].

It is well accepted that the lattice theoretic ideas are the most important technical contributions of the large literature on infinite-dimensional modern mathematical finance (for example lattice theoretic ideas in general equilibrium theory). However, ordered vector spaces that are not lattice ordered arise naturally in models of portfolio trading. Moreover, if available securities have smooth payoffs, then the portfolio space is never a vector lattice. It should be pointed out that since call and put options are vector lattice operations in the space of contingent claims, their replication by available securities requires a vector lattice structure in the portfolio space. There is a large literature on vector lattice theory related to mathematical economics; see for instance [1–7, 17–20, 22, 23].

On the other hand, there is an obvious need for properly structured high performance computational methods on vector lattices. Moreover, the main concern is to describe, in

©2012 Katsikis, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. ©2012 Katsikis, licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 2 180 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 3 MATLAB Aided Option Replication <sup>3</sup>

computational terms, and then solve problems arising from mathematical economics such as portfolio insurance and option replication. A lot of work in this area has been done in [11, 12, 14**?** –16].

In this chapter, all the numerical tasks have been performed using the Matlab 7.8 (R2009a) environment on an Intel(R) Pentium(R) Dual CPU T2310 @ 1.46 GHz 1.47 GHz 32-bit system with 2 GB of RAM memory running on the Windows Vista Home Premium Operating System.

**2. Basic results in the theory of positive bases and projection bases of R***<sup>m</sup>* In this section, a brief introduction is provided to the theory of vector lattices of **R***m*. Furthermore, we present some basic results related to the theory of positive bases and

Initially, we recall some definitions and notation from the vector lattice theory. Let **R***<sup>m</sup>* = {*<sup>x</sup>* = (*x*(1), *<sup>x</sup>*(2), ..., *<sup>x</sup>*(*m*))|*x*(*i*) <sup>∈</sup> **<sup>R</sup>**, *for each i*}, where we view **<sup>R</sup>***<sup>m</sup>* as an ordered space. The

*x* ≤ *y if and only if x*(*i*) ≤ *y*(*i*), *for each i* = 1, ..., *m*.

that *X* is a vector subspace of **R***<sup>m</sup>* then *X* ordered by the pointwise ordering is an *ordered*

basis of **<sup>R</sup>***m*. A point *<sup>x</sup>* <sup>∈</sup> **<sup>R</sup>***<sup>m</sup>* is an *upper bound, (resp. lower bound)* of a subset *<sup>S</sup>* <sup>⊆</sup> **<sup>R</sup>***<sup>m</sup>* if and only if *y* ≤ *x* (resp. *x* ≤ *y*), for all *y* ∈ *S*. For a two-point set *S* = {*x*, *y*}, we denote by *x* ∨ *y* (resp. *x* ∧ *y*) the *supremum* of *S* i.e., its least upper bound (resp. the *infimum* of *S* i.e., its greatest lower bound). Thus, *x* ∨ *y* (resp. *x* ∧ *y*) is the componentwise maximum (resp.

(*x* ∨ *y*)(*i*) = max{*x*(*i*), *y*(*i*)} ((*x* ∧ *y*)(*i*) = min{*x*(*i*), *y*(*i*)}), *for all i* = 1, ..., *m*.

For any *<sup>x</sup>* = (*x*(1), *<sup>x</sup>*(2), ..., *<sup>x</sup>*(*m*)) <sup>∈</sup> **<sup>R</sup>***m*, the set *supp*(*x*) = {*i*|*x*(*i*) �<sup>=</sup> <sup>0</sup>} is the *support* of *<sup>x</sup>*. The

An ordered subspace *<sup>Z</sup>* of **<sup>R</sup>***<sup>m</sup>* is a *vector sublattice* or a *Riesz subspace* of **<sup>R</sup>***<sup>m</sup>* if for any *<sup>x</sup>*, *<sup>y</sup>* <sup>∈</sup> *<sup>Z</sup>*

Assume that *<sup>X</sup>* is an ordered subspace of **<sup>R</sup>***<sup>m</sup>* and *<sup>B</sup>* <sup>=</sup> {*b*1, *<sup>b</sup>*2, ..., *<sup>b</sup>μ*} is a basis for *<sup>X</sup>*. Then *<sup>B</sup>* is a *positive basis* of *X* if for each *x* ∈ *X* it holds: *x* is positive if and only if its coefficients in the basis *B* are positive. In other words, *B* is a positive basis of *X* if the positive cone *X*+ of *X* has

*i* = 1, 2, ..., *μ*. A positive basis *B* = {*b*1, *b*2, ..., *bμ*} is a *partition of the unit* if the vectors *bi* have

Recall that a nonzero element *x*<sup>0</sup> of *X*<sup>+</sup> is an extremal point of *X*<sup>+</sup> if, for any *x* ∈ *X*, 0 ≤ *x* ≤ *x*<sup>0</sup> implies *x* = *λx*0, for a real number *λ*. Since each element *bi* of the positive basis of *X* is an extremal point of *X*+, a positive basis of *X* is unique in the sense of positive multiples.

*λibi*|*λ<sup>i</sup>* ≥ 0, *for each i*}.

*μ* ∑ *i*=1

<sup>+</sup> <sup>=</sup> {*<sup>x</sup>* <sup>∈</sup> **<sup>R</sup>***m*|*x*(*i*) <sup>≥</sup> 0, *for each i*} and if we suppose

*<sup>i</sup>*=<sup>1</sup> *�ibi* we have *x* ≤ *y* if and only if *λ<sup>i</sup>* ≤ *�<sup>i</sup>* for each

<sup>+</sup>. By {*e*1,*e*2, ...,*em*} we shall denote the usual

MATLAB Aided Option Replication 181

projection bases of **R***m*.

**2.1. Preliminaries and notation**

*pointwise order* relation in **R***<sup>m</sup>* is defined by

The positive cone of **R***<sup>m</sup>* is defined by **R***<sup>m</sup>*

minimum) of *x* and *y* defined by

the form,

Then, for any *<sup>x</sup>* = <sup>∑</sup>*<sup>μ</sup>*

disjoint supports and <sup>∑</sup>*<sup>μ</sup>*

*subspace* of **<sup>R</sup>***m*, with positive cone *<sup>X</sup>*<sup>+</sup> <sup>=</sup> *<sup>X</sup>* <sup>∩</sup> **<sup>R</sup>***<sup>m</sup>*

vectors *<sup>x</sup>*, *<sup>y</sup>* <sup>∈</sup> **<sup>R</sup>***<sup>m</sup>* have *disjoint supports* if *supp*(*x*) <sup>∩</sup> *supp*(*y*) = <sup>∅</sup>.

the supremum and the infimum of the set {*x*, *<sup>y</sup>*} in **<sup>R</sup>***<sup>m</sup>* belong to *<sup>Z</sup>*.

*X*<sup>+</sup> = {*x* =

*<sup>i</sup>*=<sup>1</sup> *<sup>λ</sup>ibi* and *<sup>y</sup>* <sup>=</sup> <sup>∑</sup>*<sup>μ</sup>*

*<sup>i</sup>*=<sup>1</sup> *bi* = 1.

In this chapter, we focus to the option replication problem. We consider an incomplete market of primitive securities, meaning that some call and put options need not be marketed and our objective is to describe an efficient method for computing maximal submarkets that replicate any option. Even though, there are several important results on option replication they cannot provide a method for the determination of the replicated options. By using the theory of vector-lattices and positive bases it is provided a procedure in order to determine the set of securities with replicated options. In particular, it is shown that the union of all maximal replicated submarkets (i.e., submarkets *Y*, such that any option written on the elements of *Y* can be replicated and *Y* is as large as possible) defines a set of elements such that any option written on these elements is replicated.

In [11–14, 16], it was shown that it is possible to construct computational methods in order to efficiently compute vector sublattices and lattice-subspaces of **R***<sup>m</sup>* as well as in the general case of *C*[*a*, *b*]. In addition, these methods has been successfully applied in portfolio insurance and in completion of security markets.

Here we consider a two-period security market *X* with a finite number *m* of states and a finite number of primitive securities with payoffs in **R***<sup>m</sup>* and we construct computational methods in order to determine maximal replicated submarkets of *X* by using the theory of vector sublattices and lattice-subspaces. Moreover, in the theory of security markets it is a usual practice to take call and put options with respect to the riskless bond **1** = (1, 1, ..., 1). Then, the completion *F***1**(*X*) of *X* by options is the subspace of **R***<sup>m</sup>* generated by all options written on the elements of *<sup>X</sup>* ∪ {**1**}. Since the payoff space is **<sup>R</sup>***m*, which is a vector lattice, in the case where **1** ∈ *X* then *F***1**(*X*) is exactly the vector sublattice generated by *X*. If, in addition, *X* is a vector sublattice of **R***<sup>m</sup>* then *F***1**(*X*) = *X* therefore any option is replicated, unfortunately this situation is extremely rare.

A recent article, [15], provided a computationally efficient method for computing maximal submarkets that replicate any option. In particular, by using computational methods and techniques from [11–14] in order to determine vector sublattices and their positive bases, it is presented a procedure in order to calculate the set of securities with replicated options. The aforementioned article emphasizes the most important interrelationship between the theory of vector lattices, positive bases, projection bases and the problem of option replication.

The material in this chapter is spread out in 5 sections. Section 2 is divided in two subsections; in the first the fundamental properties of lattice-subspaces and vector sublattices are presented, whereas in the second we introduce the basic results for vector sublattices, positive bases and projection bases of **R***<sup>k</sup>* together with the solution to the problem of whether a finite collection of linearly independent, positive vectors of **R***<sup>k</sup>* generates a lattice-subspace or a vector sublattice. In section 3, there are three subsections where it is discussed the theoretical background for option replication. Also, section 3 emphasis the most important interrelationship between positive bases, projection bases and the problem of option replication. Section 4 presents an algorithm for calculating maximal submarkets that replicate any option followed by the corresponding computational approach. Conclusions and research directions are provided in Section 5.

In this chapter, all the numerical tasks have been performed using the Matlab 7.8 (R2009a) environment on an Intel(R) Pentium(R) Dual CPU T2310 @ 1.46 GHz 1.47 GHz 32-bit system with 2 GB of RAM memory running on the Windows Vista Home Premium Operating System.

## **2. Basic results in the theory of positive bases and projection bases of R***<sup>m</sup>*

In this section, a brief introduction is provided to the theory of vector lattices of **R***m*. Furthermore, we present some basic results related to the theory of positive bases and projection bases of **R***m*.

#### **2.1. Preliminaries and notation**

2

[11, 12, 14**?** –16].

written on these elements is replicated.

and in completion of security markets.

this situation is extremely rare.

and research directions are provided in Section 5.

computational terms, and then solve problems arising from mathematical economics such as portfolio insurance and option replication. A lot of work in this area has been done in

In this chapter, we focus to the option replication problem. We consider an incomplete market of primitive securities, meaning that some call and put options need not be marketed and our objective is to describe an efficient method for computing maximal submarkets that replicate any option. Even though, there are several important results on option replication they cannot provide a method for the determination of the replicated options. By using the theory of vector-lattices and positive bases it is provided a procedure in order to determine the set of securities with replicated options. In particular, it is shown that the union of all maximal replicated submarkets (i.e., submarkets *Y*, such that any option written on the elements of *Y* can be replicated and *Y* is as large as possible) defines a set of elements such that any option

In [11–14, 16], it was shown that it is possible to construct computational methods in order to efficiently compute vector sublattices and lattice-subspaces of **R***<sup>m</sup>* as well as in the general case of *C*[*a*, *b*]. In addition, these methods has been successfully applied in portfolio insurance

Here we consider a two-period security market *X* with a finite number *m* of states and a finite number of primitive securities with payoffs in **R***<sup>m</sup>* and we construct computational methods in order to determine maximal replicated submarkets of *X* by using the theory of vector sublattices and lattice-subspaces. Moreover, in the theory of security markets it is a usual practice to take call and put options with respect to the riskless bond **1** = (1, 1, ..., 1). Then, the completion *F***1**(*X*) of *X* by options is the subspace of **R***<sup>m</sup>* generated by all options written on the elements of *<sup>X</sup>* ∪ {**1**}. Since the payoff space is **<sup>R</sup>***m*, which is a vector lattice, in the case where **1** ∈ *X* then *F***1**(*X*) is exactly the vector sublattice generated by *X*. If, in addition, *X* is a vector sublattice of **R***<sup>m</sup>* then *F***1**(*X*) = *X* therefore any option is replicated, unfortunately

A recent article, [15], provided a computationally efficient method for computing maximal submarkets that replicate any option. In particular, by using computational methods and techniques from [11–14] in order to determine vector sublattices and their positive bases, it is presented a procedure in order to calculate the set of securities with replicated options. The aforementioned article emphasizes the most important interrelationship between the theory of vector lattices, positive bases, projection bases and the problem of option replication.

The material in this chapter is spread out in 5 sections. Section 2 is divided in two subsections; in the first the fundamental properties of lattice-subspaces and vector sublattices are presented, whereas in the second we introduce the basic results for vector sublattices, positive bases and projection bases of **R***<sup>k</sup>* together with the solution to the problem of whether a finite collection of linearly independent, positive vectors of **R***<sup>k</sup>* generates a lattice-subspace or a vector sublattice. In section 3, there are three subsections where it is discussed the theoretical background for option replication. Also, section 3 emphasis the most important interrelationship between positive bases, projection bases and the problem of option replication. Section 4 presents an algorithm for calculating maximal submarkets that replicate any option followed by the corresponding computational approach. Conclusions Initially, we recall some definitions and notation from the vector lattice theory. Let **R***<sup>m</sup>* = {*<sup>x</sup>* = (*x*(1), *<sup>x</sup>*(2), ..., *<sup>x</sup>*(*m*))|*x*(*i*) <sup>∈</sup> **<sup>R</sup>**, *for each i*}, where we view **<sup>R</sup>***<sup>m</sup>* as an ordered space. The *pointwise order* relation in **R***<sup>m</sup>* is defined by

$$\text{If } \mathbf{x} \le y \text{ if and only if } \mathbf{x}(\mathbf{i}) \le y(\mathbf{i}) \text{, for each } \mathbf{i} = 1, \dots, m.$$

The positive cone of **R***<sup>m</sup>* is defined by **R***<sup>m</sup>* <sup>+</sup> <sup>=</sup> {*<sup>x</sup>* <sup>∈</sup> **<sup>R</sup>***m*|*x*(*i*) <sup>≥</sup> 0, *for each i*} and if we suppose that *X* is a vector subspace of **R***<sup>m</sup>* then *X* ordered by the pointwise ordering is an *ordered subspace* of **<sup>R</sup>***m*, with positive cone *<sup>X</sup>*<sup>+</sup> <sup>=</sup> *<sup>X</sup>* <sup>∩</sup> **<sup>R</sup>***<sup>m</sup>* <sup>+</sup>. By {*e*1,*e*2, ...,*em*} we shall denote the usual basis of **<sup>R</sup>***m*. A point *<sup>x</sup>* <sup>∈</sup> **<sup>R</sup>***<sup>m</sup>* is an *upper bound, (resp. lower bound)* of a subset *<sup>S</sup>* <sup>⊆</sup> **<sup>R</sup>***<sup>m</sup>* if and only if *y* ≤ *x* (resp. *x* ≤ *y*), for all *y* ∈ *S*. For a two-point set *S* = {*x*, *y*}, we denote by *x* ∨ *y* (resp. *x* ∧ *y*) the *supremum* of *S* i.e., its least upper bound (resp. the *infimum* of *S* i.e., its greatest lower bound). Thus, *x* ∨ *y* (resp. *x* ∧ *y*) is the componentwise maximum (resp. minimum) of *x* and *y* defined by

$$(\mathbf{x} \vee \mathbf{y})(i) = \max\{\mathbf{x}(i), \mathbf{y}(i)\} \left( (\mathbf{x} \wedge \mathbf{y})(i) = \min\{\mathbf{x}(i), \mathbf{y}(i)\} \right), \text{ for all } i = 1, \dots, m.$$

For any *<sup>x</sup>* = (*x*(1), *<sup>x</sup>*(2), ..., *<sup>x</sup>*(*m*)) <sup>∈</sup> **<sup>R</sup>***m*, the set *supp*(*x*) = {*i*|*x*(*i*) �<sup>=</sup> <sup>0</sup>} is the *support* of *<sup>x</sup>*. The vectors *<sup>x</sup>*, *<sup>y</sup>* <sup>∈</sup> **<sup>R</sup>***<sup>m</sup>* have *disjoint supports* if *supp*(*x*) <sup>∩</sup> *supp*(*y*) = <sup>∅</sup>.

An ordered subspace *<sup>Z</sup>* of **<sup>R</sup>***<sup>m</sup>* is a *vector sublattice* or a *Riesz subspace* of **<sup>R</sup>***<sup>m</sup>* if for any *<sup>x</sup>*, *<sup>y</sup>* <sup>∈</sup> *<sup>Z</sup>* the supremum and the infimum of the set {*x*, *<sup>y</sup>*} in **<sup>R</sup>***<sup>m</sup>* belong to *<sup>Z</sup>*.

Assume that *<sup>X</sup>* is an ordered subspace of **<sup>R</sup>***<sup>m</sup>* and *<sup>B</sup>* <sup>=</sup> {*b*1, *<sup>b</sup>*2, ..., *<sup>b</sup>μ*} is a basis for *<sup>X</sup>*. Then *<sup>B</sup>* is a *positive basis* of *X* if for each *x* ∈ *X* it holds: *x* is positive if and only if its coefficients in the basis *B* are positive. In other words, *B* is a positive basis of *X* if the positive cone *X*+ of *X* has the form,

$$X\_+ = \{ \mathbf{x} = \sum\_{i=1}^{\mu} \lambda\_i b\_i | \lambda\_i \ge 0, \text{ for each } i \}.$$

Then, for any *<sup>x</sup>* = <sup>∑</sup>*<sup>μ</sup> <sup>i</sup>*=<sup>1</sup> *<sup>λ</sup>ibi* and *<sup>y</sup>* <sup>=</sup> <sup>∑</sup>*<sup>μ</sup> <sup>i</sup>*=<sup>1</sup> *�ibi* we have *x* ≤ *y* if and only if *λ<sup>i</sup>* ≤ *�<sup>i</sup>* for each *i* = 1, 2, ..., *μ*. A positive basis *B* = {*b*1, *b*2, ..., *bμ*} is a *partition of the unit* if the vectors *bi* have disjoint supports and <sup>∑</sup>*<sup>μ</sup> <sup>i</sup>*=<sup>1</sup> *bi* = 1.

Recall that a nonzero element *x*<sup>0</sup> of *X*<sup>+</sup> is an extremal point of *X*<sup>+</sup> if, for any *x* ∈ *X*, 0 ≤ *x* ≤ *x*<sup>0</sup> implies *x* = *λx*0, for a real number *λ*. Since each element *bi* of the positive basis of *X* is an extremal point of *X*+, a positive basis of *X* is unique in the sense of positive multiples.

#### 4 182 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 3 MATLAB Aided Option Replication <sup>5</sup>

The existence of positive bases is not always ensured, but in the case where *X* is a vector sublattice of **<sup>R</sup>***<sup>m</sup>* then *<sup>X</sup>* always has a positive basis. If *<sup>B</sup>* <sup>=</sup> {*b*1, *<sup>b</sup>*2, ..., *<sup>b</sup>μ*} is a positive basis for a vector sublattice *X* the lattice operations in *X*, namely *x* ∨ *y* for the supremum and *x* ∧ *y* for the infimum of the set {*x*, *y*} in *X*, are given by

(*ii*) *Let μ* > *n. If Is* = *β*−1(*Ps*)*, and*

<sup>1</sup>, *P*� <sup>2</sup>, ..., *P*�

*then*

*R*(*γ*) = {*P*�

*(i) Z* = *L* ⊕ [*zr*+1, ..., *zμ*],

 *<sup>i</sup>* + *b*� *i , with b* 

Then, by Theorem 2, if

*(iii) If bi* = *b*

{*b* 1, *b* 2, ..., *b*

*(ii)* {*br*+1, *br*+2, ..., *bμ*} = {2*zr*+1, 2*zr*+2, ..., 2*zμ*},

*basis for L which is given by the formula*

*coordinates of x in the projection basis, i.e.,*

*<sup>i</sup>* ∈ *L and b*�

(*b* 1, *b* 2, ..., *b*

> *x* = *μ* ∑ *i*=1

positive vectors of **<sup>R</sup>***<sup>m</sup>* and {*d*1, *<sup>d</sup>*2, ..., *<sup>d</sup>μ*} is a positive basis for *<sup>Z</sup>*.

(*d* 1, *d* 2, ..., *d*

*xs* = ∑ *i*∈*Is*

*is the vector sublattice generated by x*1, *x*2, ..., *xn and* dim *Z* = *μ*.

*<sup>μ</sup>*}*. Then, the relation*

*where B is the μ* × *μ matrix with columns the vectors P*�

�*h*(*i*)�1*ei*, *s* = *n* + 1, *n* + 2, ..., *μ*,

(*b*1, *b*2, ..., *bμ*)*<sup>T</sup>* = *B*−1(*x*1, *x*2, ..., *xμ*)*<sup>T</sup>* (2)

*<sup>μ</sup>, defines a positive basis for Z.*

MATLAB Aided Option Replication 183

 1, *b* 2, ..., *b*

*<sup>r</sup>*} *is a*

<sup>1</sup>, *P*� <sup>2</sup>, ..., *P*�

*<sup>i</sup>* ∈ [*zr*+1, ..., *zμ*]*, for each i* = 1, 2, ...,*r, then* {*b*

*r* ∑ *i*=1 *λib i*

*<sup>r</sup>*)*<sup>T</sup>* <sup>=</sup> *<sup>A</sup>*−1(*z*1, *<sup>z</sup>*2, ..., *zr*)*T*,

*<sup>r</sup>*} *is called the projection basis of L and has the property: The r first coordinates of*

*where A is the r* × *r matrix whose ith column is the vector Pi, for i* = 1, 2, ...,*r. This basis,*

*any element x* ∈ *L expressed in terms of the basis* {*b*1, *b*2, ..., *bμ*} *coincide with the corresponding*

*λibi* ∈ *L* ⇒ *x* =

*<sup>r</sup>*)*<sup>T</sup>* <sup>=</sup> *<sup>A</sup>*−1(*z*1, *<sup>z</sup>*2, ..., *zr*)*T*,

Suppose that *Z* is the sublattice generated by a collection *z*1, *z*2, ..., *zr* of linearly independent,

*Z* = [*x*1, *x*2, ..., *xn*, *xn*+1, ..., *xμ*],

*For a positive basis* {*b*1, *b*2, ..., *bμ*} *of Z, consider the basic function γ of* {*x*1, *x*2, ..., *xμ*} *with range,*

The notion of the projection basis is important for our study. Furthermore, in the following, we are interested for a projection basis that corresponds to a positive basis. Let {*z*1, *z*2, ..., *zr*} be a set of linearly independent and positive vectors of **R***<sup>m</sup>* then by using Theorem 1 we construct the sublattice *Z* of **R***<sup>m</sup>* generated by these vectors. If dim(*Z*) = *μ*, by Theorem 1, a positive basis {*b*1, *b*2, ..., *bμ*} of *Z* can be determined. The basic result for calculating the projection basis

**Theorem 2.** *Suppose that β is the basic function of the vectors* {*z*1, *z*2, ..., *zr*} *and P*1, *P*2, ..., *Pr*, *Pr*+1, ..., *P<sup>μ</sup> is an enumeration of the range of β such that the first r vectors P*1, *P*2, ..., *Pr are linearly independent and suppose also that zr*+1, ..., *z<sup>μ</sup> are the new vectors constructed in Theorem 1. If*

that corresponds to the positive basis {*b*1, *b*2, ..., *bμ*} of *Z* is the following theorem.

*L* = [*z*1, *z*2, ..., *zr*] *is the subspace of* **R***<sup>m</sup> generated by the vectors z*1, *z*2, ..., *zr then,*

$$\begin{aligned} \text{in } \boldsymbol{x} \lor \boldsymbol{y} &= \sum\_{i=1}^{\mu} \max\{\lambda\_{i\prime}, \varrho\_{i}\} b\_{i} \text{ and } \boldsymbol{x} \land \boldsymbol{y} = \sum\_{i=1}^{\mu} \min\{\lambda\_{i\prime}, \varrho\_{i}\} b\_{i\prime}, \\\ \boldsymbol{x} \land \boldsymbol{y} &= \sum\_{i=1}^{\mu} \max\{\lambda\_{i\prime}, \varrho\_{i}\} b\_{i\prime} \end{aligned}$$

for each *<sup>x</sup>* = <sup>∑</sup>*<sup>μ</sup> <sup>i</sup>*=<sup>1</sup> *<sup>λ</sup>ibi*, *<sup>y</sup>* <sup>=</sup> <sup>∑</sup>*<sup>μ</sup> <sup>i</sup>*=<sup>1</sup> *�ibi* ∈ *X*.

Suppose that *L* is a finite dimensional subspace of *C*(Ω) generated by a set {*z*1, *z*2, ..., *zr*} of linearly independent positive vectors of *C*(Ω). If *Z* is the sublattice of *C*(Ω) generated by *L* and {*b*1, ..., *<sup>b</sup>μ*} is a positive basis for *<sup>Z</sup>* (*<sup>μ</sup>* <sup>=</sup> *dim*(*Z*)) then, a *projection basis* { ˜*b*1, ˜*b*2, ..., ˜ *br*} of *Z* is a basis for *L* such that its elements are projections of the elements of the positive basis {*b*1, ..., *<sup>b</sup>μ*}. In our current work we consider that <sup>Ω</sup> <sup>=</sup> {1, 2, ..., *<sup>m</sup>*} hence *<sup>C</sup>*(Ω) = **<sup>R</sup>***m*.

For an extensive presentation of vector sublattices as well as for notation not defined here we refer to [11–13, 15**?** , 16] and the references therein.

#### **2.2. Theoretical background**

In this section we present theoretical results for positive bases and projection bases in **R***m*. Given a collection *x*1, *x*2, ..., *xn* of linearly independent, positive vectors of **R***<sup>m</sup>* we define the function,

$$h: \{1, 2, \ldots, m\} \to \mathbb{R}^n \text{ such that } h(i) = (\mathfrak{x}\_1(i), \mathfrak{x}\_2(i), \ldots, \mathfrak{x}\_n(i))$$

and the function,

$$\beta: \{1, 2, \ldots, m\} \to \mathbb{R}^n \text{ such that } \beta(i) = \frac{h(i)}{||h(i)||\_1} \tag{1}$$

for each *i* ∈ {1, 2, ..., *m*} with �*h*(*i*)�<sup>1</sup> �= 0. We shall refer to *β* as the *basic function* of the vectors *x*1, *x*2, ..., *xn*. The set

$$R(\beta) = \{\beta(i)|i = 1, 2, \dots, m, \text{ with } ||h(i)||\_1 \neq 0\},$$

is the *range* of the basic function and the *cardinal number*, *cardR*(*β*), of *R*(*β*) is the number of different elements of *R*(*β*).

Suppose that *Z* denotes the sublattice of **R***<sup>m</sup>* generated by *X* = [*x*1, *x*2, ..., *xn*]. We shall denote by *P*1, *P*2, ..., *Pn*, *Pn*+1, ..., *P<sup>μ</sup>* an enumeration of *R*(*β*) such that the first *n* vertices *P*1, *P*2, ..., *Pn* are linearly independent and *μ* = dim(*Z*). Note that such an enumeration always exists. The notation, *A<sup>T</sup>* stands for the transpose of a matrix *A*. A procedure in order to construct the sublattice *Z* is given by the following theorem.

**Theorem 1.** *Suppose that the above assumptions are satisfied. Then,*

(*i*) *X is a vector sublattice of* **R***<sup>m</sup> if and only if R*(*β*) *has exactly n points (i.e., μ* = *n). Then a positive basis b*1, *b*2, ..., *bn for X is defined by the formula*

$$(\boldsymbol{b}\_1, \boldsymbol{b}\_2, \dots, \boldsymbol{b}\_n)^T = A^{-1} (\boldsymbol{\chi}\_1, \boldsymbol{\chi}\_2, \dots, \boldsymbol{\chi}\_n)^T$$

*where A is the n* × *n matrix whose ith column is the vector Pi*, *for each i* = 1, 2, ..., *n*. *It is clear that in such a case Z and X coincide.*

(*ii*) *Let μ* > *n. If Is* = *β*−1(*Ps*)*, and*

$$\mathfrak{x}\_{\mathbf{s}} = \sum\_{i \in I\_{\mathbf{s}}} ||h(i)||\_1 e\_{i\prime} \text{ s} = n + 1, n + 2, \dots, \mu\_{\mathbf{s}}$$

*then*

4

for each *<sup>x</sup>* = <sup>∑</sup>*<sup>μ</sup>*

function,

and the function,

*x*1, *x*2, ..., *xn*. The set

different elements of *R*(*β*).

sublattice *Z* is given by the following theorem.

*basis b*1, *b*2, ..., *bn for X is defined by the formula*

*that in such a case Z and X coincide.*

**Theorem 1.** *Suppose that the above assumptions are satisfied. Then,*

The existence of positive bases is not always ensured, but in the case where *X* is a vector sublattice of **<sup>R</sup>***<sup>m</sup>* then *<sup>X</sup>* always has a positive basis. If *<sup>B</sup>* <sup>=</sup> {*b*1, *<sup>b</sup>*2, ..., *<sup>b</sup>μ*} is a positive basis for a vector sublattice *X* the lattice operations in *X*, namely *x* ∨ *y* for the supremum and *x* ∧ *y* for

> *μ* ∑ *i*=1

min{*λi*, *�i*}*bi*,

�*h*(*i*)�<sup>1</sup>

*br*} of

(1)

max{*λi*, *�i*}*bi and x* ∧ *y* =

Suppose that *L* is a finite dimensional subspace of *C*(Ω) generated by a set {*z*1, *z*2, ..., *zr*} of linearly independent positive vectors of *C*(Ω). If *Z* is the sublattice of *C*(Ω) generated by *L* and {*b*1, ..., *<sup>b</sup>μ*} is a positive basis for *<sup>Z</sup>* (*<sup>μ</sup>* <sup>=</sup> *dim*(*Z*)) then, a *projection basis* { ˜*b*1, ˜*b*2, ..., ˜

*Z* is a basis for *L* such that its elements are projections of the elements of the positive basis

For an extensive presentation of vector sublattices as well as for notation not defined here we

In this section we present theoretical results for positive bases and projection bases in **R***m*. Given a collection *x*1, *x*2, ..., *xn* of linearly independent, positive vectors of **R***<sup>m</sup>* we define the

*<sup>h</sup>* : {1, 2, ..., *<sup>m</sup>*} → **<sup>R</sup>***<sup>n</sup> such that h*(*i*)=(*x*1(*i*), *<sup>x</sup>*2(*i*), ..., *xn*(*i*))

*<sup>β</sup>* : {1, 2, ..., *<sup>m</sup>*} → **<sup>R</sup>***<sup>n</sup> such that <sup>β</sup>*(*i*) = *<sup>h</sup>*(*i*)

for each *i* ∈ {1, 2, ..., *m*} with �*h*(*i*)�<sup>1</sup> �= 0. We shall refer to *β* as the *basic function* of the vectors

*R*(*β*) = {*β*(*i*)|*i* = 1, 2, ..., *m*, *with* �*h*(*i*)�<sup>1</sup> �= 0}, is the *range* of the basic function and the *cardinal number*, *cardR*(*β*), of *R*(*β*) is the number of

Suppose that *Z* denotes the sublattice of **R***<sup>m</sup>* generated by *X* = [*x*1, *x*2, ..., *xn*]. We shall denote by *P*1, *P*2, ..., *Pn*, *Pn*+1, ..., *P<sup>μ</sup>* an enumeration of *R*(*β*) such that the first *n* vertices *P*1, *P*2, ..., *Pn* are linearly independent and *μ* = dim(*Z*). Note that such an enumeration always exists. The notation, *A<sup>T</sup>* stands for the transpose of a matrix *A*. A procedure in order to construct the

(*i*) *X is a vector sublattice of* **R***<sup>m</sup> if and only if R*(*β*) *has exactly n points (i.e., μ* = *n). Then a positive*

(*b*1, *b*2, ..., *bn*)*<sup>T</sup>* = *A*−1(*x*1, *x*2, ..., *xn*)*T*,

*where A is the n* × *n matrix whose ith column is the vector Pi*, *for each i* = 1, 2, ..., *n*. *It is clear*

{*b*1, ..., *<sup>b</sup>μ*}. In our current work we consider that <sup>Ω</sup> <sup>=</sup> {1, 2, ..., *<sup>m</sup>*} hence *<sup>C</sup>*(Ω) = **<sup>R</sup>***m*.

the infimum of the set {*x*, *y*} in *X*, are given by

*μ* ∑ *i*=1

*<sup>i</sup>*=<sup>1</sup> *�ibi* ∈ *X*.

*x* ∨ *y* =

*<sup>i</sup>*=<sup>1</sup> *<sup>λ</sup>ibi*, *<sup>y</sup>* <sup>=</sup> <sup>∑</sup>*<sup>μ</sup>*

refer to [11–13, 15**?** , 16] and the references therein.

**2.2. Theoretical background**

$$Z = [\mathfrak{x}\_1, \mathfrak{x}\_2, \dots, \mathfrak{x}\_{n'} \mathfrak{x}\_{n+1}, \dots, \mathfrak{x}\_{\mu}]\_{\mathsf{A}}$$

*is the vector sublattice generated by x*1, *x*2, ..., *xn and* dim *Z* = *μ*. *For a positive basis* {*b*1, *b*2, ..., *bμ*} *of Z, consider the basic function γ of* {*x*1, *x*2, ..., *xμ*} *with range, R*(*γ*) = {*P*� <sup>1</sup>, *P*� <sup>2</sup>, ..., *P*� *<sup>μ</sup>*}*. Then, the relation*

$$(b\_1, b\_2, \ldots, b\_{\mu})^T = B^{-1} (\mathbf{x}\_1, \mathbf{x}\_2, \ldots, \mathbf{x}\_{\mu})^T \tag{2}$$

*where B is the μ* × *μ matrix with columns the vectors P*� <sup>1</sup>, *P*� <sup>2</sup>, ..., *P*� *<sup>μ</sup>, defines a positive basis for Z.*

The notion of the projection basis is important for our study. Furthermore, in the following, we are interested for a projection basis that corresponds to a positive basis. Let {*z*1, *z*2, ..., *zr*} be a set of linearly independent and positive vectors of **R***<sup>m</sup>* then by using Theorem 1 we construct the sublattice *Z* of **R***<sup>m</sup>* generated by these vectors. If dim(*Z*) = *μ*, by Theorem 1, a positive basis {*b*1, *b*2, ..., *bμ*} of *Z* can be determined. The basic result for calculating the projection basis that corresponds to the positive basis {*b*1, *b*2, ..., *bμ*} of *Z* is the following theorem.

**Theorem 2.** *Suppose that β is the basic function of the vectors* {*z*1, *z*2, ..., *zr*} *and P*1, *P*2, ..., *Pr*, *Pr*+1, ..., *P<sup>μ</sup> is an enumeration of the range of β such that the first r vectors P*1, *P*2, ..., *Pr are linearly independent and suppose also that zr*+1, ..., *z<sup>μ</sup> are the new vectors constructed in Theorem 1. If L* = [*z*1, *z*2, ..., *zr*] *is the subspace of* **R***<sup>m</sup> generated by the vectors z*1, *z*2, ..., *zr then,*


$$(\tilde{b\_1}, \tilde{b\_2}, \dots, \tilde{b\_r})^T = A^{-1} (z\_{1'} z\_{2'} \dots z\_r)^T$$

*where A is the r* × *r matrix whose ith column is the vector Pi, for i* = 1, 2, ...,*r. This basis,* {*b* 1, *b* 2, ..., *b <sup>r</sup>*} *is called the projection basis of L and has the property: The r first coordinates of any element x* ∈ *L expressed in terms of the basis* {*b*1, *b*2, ..., *bμ*} *coincide with the corresponding coordinates of x in the projection basis, i.e.,*

$$\mathfrak{x} = \sum\_{i=1}^{\mu} \lambda\_i b\_i \in L \Rightarrow \mathfrak{x} = \sum\_{i=1}^{r} \lambda\_i \tilde{b}\_i$$

Suppose that *Z* is the sublattice generated by a collection *z*1, *z*2, ..., *zr* of linearly independent, positive vectors of **<sup>R</sup>***<sup>m</sup>* and {*d*1, *<sup>d</sup>*2, ..., *<sup>d</sup>μ*} is a positive basis for *<sup>Z</sup>*.

Then, by Theorem 2, if

$$(\tilde{d}\_1, \tilde{d}\_2, \dots, \tilde{d}\_r)^T = A^{-1} (z\_1, z\_2, \dots, z\_r)^T, \dots$$

where *A* is the *r* × *r* matrix whose *i*th column is the vector *Pi*, for *i* = 1, 2, ...,*r* then {*d* 1, *d* 2, ..., *d <sup>r</sup>*} is the projection basis of *L* = [*z*1, *z*2, ..., *zr*]. The projection basis {*d* 1, *d* 2, ..., *d r*} is called the *projection basis of X corresponding to the basis* {*d*1, *d*2, ..., *dμ*}. The following proposition allows us to determine a positive basis and its corresponding projection basis. Moreover, the calculated positive basis is a partition of the unit.

**3.1. The economic model**

portfolio.

identity

from [16].

which is called *put-call parity*.

*x* − *α***1** = *c*(*x*, *α*) − *p*(*x*, *α*).

**3.2. Completion of security markets**

*s* = 1, 2, ..., *m*, while at *t* = 0 the state is known to be *s* = 0.

*A* =

*xi*(*j*)

future payoffs of *x*1, *x*2, ..., *xn* are collected in a matrix

In our economy there are two time periods, *t* = 0, 1, where *t* = 0 denotes the present and *t* = 1 denotes the future. We consider that at *t* = 1 we have a finite number of states indexed by

Suppose that, agents trade *x*1, *x*2, ..., *xn* non-redundant securities in period *t* = 0, then the

*<sup>j</sup>*=1,2,...,*<sup>m</sup>*

where *xi*(*j*) is the payoff of one unit of security *i* in state *j*. In other words, *A* is the matrix whose columns are the non-redundant security vectors *x*1, *x*2, ..., *xn*. It is clear that the matrix *A* is of full rank and the *asset span* is denoted by *X* = *Span*(*A*). So, *X* is the vector subspace of **R***<sup>m</sup>* generated by the vectors *xi*. That is, *X* consists of those income streams that can be generated by trading on the financial market. A *portfolio* is a column vector *θ* = (*θ*1, *θ*2, ..., *θn*)*<sup>T</sup>* of **<sup>R</sup>***<sup>n</sup>* and the *payoff* of a portfolio *<sup>θ</sup>* is the vector *<sup>x</sup>* <sup>=</sup> *<sup>A</sup><sup>θ</sup>* <sup>∈</sup> **<sup>R</sup>***<sup>m</sup>* which offers payoff *<sup>x</sup>*(*i*) in state *i*, where *i* = 1, ..., *m*. A vector in **R***m*, is said to be *marketed* or *replicated* if *x* is the payoff of some portfolio *θ* (called the *replicating portfolio* of *x*), or equivalently if *x* ∈ *X*. If *m* = *n*, then markets are said to be *complete* and the asset span coincides with the space **R***m*. On the other hand, if *n* < *m*, the markets are *incomplete* and some state contingent claim cannot be replicated by a

In the following, we shall denote by **1** the *riskless bond* i.e., the vector **1** = (1, 1, ..., 1). A vector *x* is a *binary vector* if *x* �= **0** = (0, 0, ..., 0), *x* �= **1** and *x*(*i*) = 0 or *x*(*i*) = 1, for any *<sup>i</sup>*. The *call option* written on the vector *<sup>x</sup>* <sup>∈</sup> **<sup>R</sup>***<sup>m</sup>* with exercise price *<sup>α</sup>* is the vector *<sup>c</sup>*(*x*, *<sup>a</sup>*) = (*<sup>x</sup>* <sup>−</sup> *<sup>α</sup>***1**)<sup>+</sup> = (*<sup>x</sup>* <sup>−</sup> *<sup>α</sup>***1**) <sup>∨</sup> **<sup>0</sup>**. The *put option* written on the vector *<sup>x</sup>* <sup>∈</sup> **<sup>R</sup>***<sup>m</sup>* with exercise price *<sup>α</sup>* is the vector *<sup>p</sup>*(*x*, *<sup>a</sup>*)=(*α***<sup>1</sup>** <sup>−</sup> *<sup>x</sup>*)<sup>+</sup> = (*α***<sup>1</sup>** <sup>−</sup> *<sup>x</sup>*) <sup>∨</sup> **<sup>0</sup>**. If *<sup>y</sup>* is an element of a Riesz space then the following lattice identities hold, *<sup>y</sup>* <sup>=</sup> *<sup>y</sup>*<sup>+</sup> <sup>−</sup> *<sup>y</sup>*<sup>−</sup> and *<sup>y</sup>*<sup>−</sup> = (−*y*)+. It is clear that *<sup>x</sup>* <sup>−</sup> *<sup>α</sup>***<sup>1</sup>** <sup>=</sup> (*<sup>x</sup>* <sup>−</sup> *<sup>α</sup>***1**)<sup>+</sup> <sup>−</sup> (*<sup>x</sup>* <sup>−</sup> *<sup>α</sup>***1**)<sup>−</sup> = (*<sup>x</sup>* <sup>−</sup> *<sup>α</sup>***1**)<sup>+</sup> <sup>−</sup> (*α***<sup>1</sup>** <sup>−</sup> *<sup>x</sup>*)<sup>+</sup> <sup>=</sup> *<sup>c</sup>*(*x*, *<sup>α</sup>*) <sup>−</sup> *<sup>p</sup>*(*x*, *<sup>α</sup>*). Therefore we have the

*x* − *a***1** = *c*(*x*, *a*) − *p*(*x*, *a*),

If both *c*(*x*, *α*) > 0 and *p*(*x*, *α*) > 0, we say that the call option *c*(*x*, *α*) and the put option *p*(*x*, *α*) are *non trivial* and the exercise price *α* is a *non trivial exercise price* of *x*. If *c*(*x*, *α*) and *p*(*x*, *α*) belong to *X* then we say that *c*(*x*, *α*) and *p*(*x*, *α*) are *replicated*. If we suppose that **1** ∈ *X* and at least one of *c*(*x*, *a*), *p*(*x*, *a*) is replicated, then both of them are replicated since,

We shall discuss the problem of completion by options of a two-period security market in which the space of marketed securities is a subspace of **R***m*. The present study involves vector sublattices generated by a subset *B* of **R***<sup>m</sup>* of positive, linearly independent vectors. A computational solution to this problem is provided by using the SUBlat Matlab function

*<sup>i</sup>*=1,2,...*<sup>n</sup>* <sup>∈</sup> **<sup>R</sup>***m*×*<sup>n</sup>*

MATLAB Aided Option Replication 185

**Proposition 1.** *Suppose that* {*di*} *is the basis of Z given by equation (2) of Theorem 1 and* {*d <sup>i</sup>*} *is the projection basis of L* = [*z*1, *<sup>z</sup>*2, ..., *zr*] *corresponding to the basis* {*di*}. *Then* {*bi* <sup>=</sup> *di* �*di*�<sup>∞</sup> <sup>|</sup>*<sup>i</sup>* <sup>=</sup> 1, 2, ..., *<sup>μ</sup>*} *is the positive basis of Z which is a partition of the unit and* {*b <sup>i</sup>* = *<sup>d</sup> i* �*di*�<sup>∞</sup> <sup>|</sup>*<sup>i</sup>* <sup>=</sup> 1, 2, ..., *<sup>n</sup>*} *is the projection basis of L corresponding to the basis* {*bi*} *of Z.*

In the following, we shall denote by **1** the vector **1** = (1, 1, ..., 1). A vector *x* is a *binary vector* if *x* �= **0** = (0, 0, ..., 0), *x* �= **1** and *x*(*i*) = 0 or *x*(*i*) = 1, for any *i*. Let {*bi*|*i* = 1, 2, ..., *μ*} be a positive basis of *Z* which is a partition of the unit and let {*b <sup>i</sup>*|*i* = 1, 2, ..., *n*} be the projection basis of *L* corresponding to the basis {*bi*}. A partition *δ* = {*σi*|*i* = 1, 2, ..., *k*} of {1, 2, ..., *n*} is *proper* if for any *<sup>r</sup>* = 1, 2, ..., *<sup>k</sup>*, the vector *wr* = <sup>∑</sup>*i*∈*σ<sup>r</sup> <sup>b</sup> <sup>i</sup>* is a binary vector with ∑*<sup>k</sup> <sup>r</sup>*=<sup>1</sup> *wr* = **1**. If there is no proper partition of {1, 2, ..., *n*} strictly finer than *δ*, then we say that *δ* is a *maximal proper partition* of {1, 2, ..., *n*}.

**Example 1.** Let {*bi*|*i* = 1, 2, 3} be a positive basis such that the corresponding projection basis is the following

$$
\widetilde{b}\_1 = (\frac{1}{2}, 1, 0, 1, 0), \\
\widetilde{b}\_2 = (\frac{1}{2}, 0, 0, 0, 1), \\
\widetilde{b}\_3 = (0, 0, 1, 0, 0).
$$

We calculate the partitions of {1, 2, 3} which are the following:

$$\delta\_1 = \{\sigma\_1, \sigma\_2\}, \text{ where } \sigma\_1 = \{1\}, \ \sigma\_2 = \{2 \ 3\}$$

$$\delta\_2 = \{\sigma\_1, \sigma\_2\}, \text{ where } \sigma\_1 = \{2\}, \ \sigma\_2 = \{3 \ 4\}$$

$$\delta\_3 = \{\sigma\_1, \sigma\_2\}, \text{ where } \sigma\_1 = \{3\}, \ \sigma\_2 = \{1 \ 2\}$$

$$\delta\_4 = \{\sigma\_1, \sigma\_2, \sigma\_3\}, \text{ where } \sigma\_1 = \{1\}, \ \sigma\_2 = \{2\}, \ \sigma\_3 = \{3\}.$$

Then *<sup>δ</sup>*<sup>1</sup> is not a proper partition since *<sup>w</sup>*<sup>1</sup> = <sup>∑</sup>*i*∈*σ*<sup>1</sup> *bi* = *b*<sup>1</sup> and *b*<sup>1</sup> is not a binary vector. Similarly, *δ*<sup>2</sup> is not a proper partition. On the other hand *δ*<sup>3</sup> is a proper partition since *w*<sup>1</sup> = ∑*i*∈*σ*<sup>1</sup> *bi* = *<sup>b</sup>*<sup>3</sup> and *<sup>b</sup>*<sup>3</sup> is a binary vector, *<sup>w</sup>*<sup>2</sup> = <sup>∑</sup>*i*∈*σ*<sup>2</sup> *bi* = *b*<sup>1</sup> +*b*<sup>2</sup> = (1, 1, 0, 1, 1) which is a binary vector and *w*<sup>1</sup> + *w*<sup>2</sup> = (1, 1, 1, 1, 1). Notice that *δ*<sup>3</sup> is strictly finer than *δ*4, hence *δ*<sup>3</sup> is a maximal proper partition of {1, 2, 3}.

#### **3. Option replication**

In this section we shall discuss the economic model of our study. Moreover, first we discuss an inductive method for calculating the completion of security markets. So, if **1** ∈ *X* then it is possible to determine a basic set of marketed securities i.e., a set of linearly independent and positive vectors of *X* and the sublattice of **R***<sup>m</sup>* generated by a basic set of marketed securities is *F***1**(*X*). Finally, *F***1**(*X*) has a positive basis which is a partition of the unit. Under these observations we present an algorithmic procedure in order to determine maximal submarkets that replicate any option.

#### **3.1. The economic model**

6

{*d* 1, *d* 2, ..., *d*

where *A* is the *r* × *r* matrix whose *i*th column is the vector *Pi*, for *i* = 1, 2, ...,*r* then

*<sup>r</sup>*} is the projection basis of *L* = [*z*1, *z*2, ..., *zr*]. The projection basis {*d*

is called the *projection basis of X corresponding to the basis* {*d*1, *d*2, ..., *dμ*}. The following proposition allows us to determine a positive basis and its corresponding projection basis.

In the following, we shall denote by **1** the vector **1** = (1, 1, ..., 1). A vector *x* is a *binary vector* if *x* �= **0** = (0, 0, ..., 0), *x* �= **1** and *x*(*i*) = 0 or *x*(*i*) = 1, for any *i*. Let {*bi*|*i* = 1, 2, ..., *μ*} be a

basis of *L* corresponding to the basis {*bi*}. A partition *δ* = {*σi*|*i* = 1, 2, ..., *k*} of {1, 2, ..., *n*} is

there is no proper partition of {1, 2, ..., *n*} strictly finer than *δ*, then we say that *δ* is a *maximal*

**Example 1.** Let {*bi*|*i* = 1, 2, 3} be a positive basis such that the corresponding projection basis

2

*δ*<sup>1</sup> = {*σ*1, *σ*2}, *where σ*<sup>1</sup> = {1}, *σ*<sup>2</sup> = {2 3}

*δ*<sup>2</sup> = {*σ*1, *σ*2}, *where σ*<sup>1</sup> = {2}, *σ*<sup>2</sup> = {3 4} *δ*<sup>3</sup> = {*σ*1, *σ*2}, *where σ*<sup>1</sup> = {3}, *σ*<sup>2</sup> = {1 2} *δ*<sup>4</sup> = {*σ*1, *σ*2, *σ*3}, *where σ*<sup>1</sup> = {1}, *σ*<sup>2</sup> = {2}, *σ*<sup>3</sup> = {3}.

Similarly, *δ*<sup>2</sup> is not a proper partition. On the other hand *δ*<sup>3</sup> is a proper partition since *w*<sup>1</sup> =

vector and *w*<sup>1</sup> + *w*<sup>2</sup> = (1, 1, 1, 1, 1). Notice that *δ*<sup>3</sup> is strictly finer than *δ*4, hence *δ*<sup>3</sup> is a maximal

In this section we shall discuss the economic model of our study. Moreover, first we discuss an inductive method for calculating the completion of security markets. So, if **1** ∈ *X* then it is possible to determine a basic set of marketed securities i.e., a set of linearly independent and positive vectors of *X* and the sublattice of **R***<sup>m</sup>* generated by a basic set of marketed securities is *F***1**(*X*). Finally, *F***1**(*X*) has a positive basis which is a partition of the unit. Under these observations we present an algorithmic procedure in order to determine maximal submarkets

 *<sup>i</sup>* = *<sup>d</sup> i*

, 0, 0, 0, 1), *b*<sup>3</sup> = (0, 0, 1, 0, 0).

*<sup>i</sup>* is a binary vector with ∑*<sup>k</sup>*

**Proposition 1.** *Suppose that* {*di*} *is the basis of Z given by equation (2) of Theorem 1 and* {*d*

*projection basis of L* = [*z*1, *<sup>z</sup>*2, ..., *zr*] *corresponding to the basis* {*di*}. *Then* {*bi* <sup>=</sup> *di*

Moreover, the calculated positive basis is a partition of the unit.

*is the positive basis of Z which is a partition of the unit and* {*b*

positive basis of *Z* which is a partition of the unit and let {*b*

, 1, 0, 1, 0), *<sup>b</sup>*<sup>2</sup> = ( <sup>1</sup>

We calculate the partitions of {1, 2, 3} which are the following:

Then *<sup>δ</sup>*<sup>1</sup> is not a proper partition since *<sup>w</sup>*<sup>1</sup> = <sup>∑</sup>*i*∈*σ*<sup>1</sup>

*bi* = *<sup>b</sup>*<sup>3</sup> and *<sup>b</sup>*<sup>3</sup> is a binary vector, *<sup>w</sup>*<sup>2</sup> = <sup>∑</sup>*i*∈*σ*<sup>2</sup>

*proper* if for any *<sup>r</sup>* = 1, 2, ..., *<sup>k</sup>*, the vector *wr* = <sup>∑</sup>*i*∈*σ<sup>r</sup> <sup>b</sup>*

*<sup>b</sup>*<sup>1</sup> = ( <sup>1</sup> 2

*basis of L corresponding to the basis* {*bi*} *of Z.*

*proper partition* of {1, 2, ..., *n*}.

proper partition of {1, 2, 3}.

**3. Option replication**

that replicate any option.

is the following

∑*i*∈*σ*<sup>1</sup>

 1, *d* 2, ..., *d r*}

> *<sup>i</sup>*} *is the*

*<sup>r</sup>*=<sup>1</sup> *wr* = **1**. If

�*di*�<sup>∞</sup> <sup>|</sup>*<sup>i</sup>* <sup>=</sup> 1, 2, ..., *<sup>μ</sup>*}

�*di*�<sup>∞</sup> <sup>|</sup>*<sup>i</sup>* <sup>=</sup> 1, 2, ..., *<sup>n</sup>*} *is the projection*

*<sup>i</sup>*|*i* = 1, 2, ..., *n*} be the projection

*bi* = *b*<sup>1</sup> and *b*<sup>1</sup> is not a binary vector.

*bi* = *b*<sup>1</sup> +*b*<sup>2</sup> = (1, 1, 0, 1, 1) which is a binary

In our economy there are two time periods, *t* = 0, 1, where *t* = 0 denotes the present and *t* = 1 denotes the future. We consider that at *t* = 1 we have a finite number of states indexed by *s* = 1, 2, ..., *m*, while at *t* = 0 the state is known to be *s* = 0.

Suppose that, agents trade *x*1, *x*2, ..., *xn* non-redundant securities in period *t* = 0, then the future payoffs of *x*1, *x*2, ..., *xn* are collected in a matrix

$$A = \left[ \mathfrak{x}\_{i}(j) \right]\_{i=1,2,\ldots,m}^{j=1,2,\ldots,m} \in \mathbb{R}^{m \times m}$$

where *xi*(*j*) is the payoff of one unit of security *i* in state *j*. In other words, *A* is the matrix whose columns are the non-redundant security vectors *x*1, *x*2, ..., *xn*. It is clear that the matrix *A* is of full rank and the *asset span* is denoted by *X* = *Span*(*A*). So, *X* is the vector subspace of **R***<sup>m</sup>* generated by the vectors *xi*. That is, *X* consists of those income streams that can be generated by trading on the financial market. A *portfolio* is a column vector *θ* = (*θ*1, *θ*2, ..., *θn*)*<sup>T</sup>* of **<sup>R</sup>***<sup>n</sup>* and the *payoff* of a portfolio *<sup>θ</sup>* is the vector *<sup>x</sup>* <sup>=</sup> *<sup>A</sup><sup>θ</sup>* <sup>∈</sup> **<sup>R</sup>***<sup>m</sup>* which offers payoff *<sup>x</sup>*(*i*) in state *i*, where *i* = 1, ..., *m*. A vector in **R***m*, is said to be *marketed* or *replicated* if *x* is the payoff of some portfolio *θ* (called the *replicating portfolio* of *x*), or equivalently if *x* ∈ *X*. If *m* = *n*, then markets are said to be *complete* and the asset span coincides with the space **R***m*. On the other hand, if *n* < *m*, the markets are *incomplete* and some state contingent claim cannot be replicated by a portfolio.

In the following, we shall denote by **1** the *riskless bond* i.e., the vector **1** = (1, 1, ..., 1). A vector *x* is a *binary vector* if *x* �= **0** = (0, 0, ..., 0), *x* �= **1** and *x*(*i*) = 0 or *x*(*i*) = 1, for any *<sup>i</sup>*. The *call option* written on the vector *<sup>x</sup>* <sup>∈</sup> **<sup>R</sup>***<sup>m</sup>* with exercise price *<sup>α</sup>* is the vector *<sup>c</sup>*(*x*, *<sup>a</sup>*) = (*<sup>x</sup>* <sup>−</sup> *<sup>α</sup>***1**)<sup>+</sup> = (*<sup>x</sup>* <sup>−</sup> *<sup>α</sup>***1**) <sup>∨</sup> **<sup>0</sup>**. The *put option* written on the vector *<sup>x</sup>* <sup>∈</sup> **<sup>R</sup>***<sup>m</sup>* with exercise price *<sup>α</sup>* is the vector *<sup>p</sup>*(*x*, *<sup>a</sup>*)=(*α***<sup>1</sup>** <sup>−</sup> *<sup>x</sup>*)<sup>+</sup> = (*α***<sup>1</sup>** <sup>−</sup> *<sup>x</sup>*) <sup>∨</sup> **<sup>0</sup>**. If *<sup>y</sup>* is an element of a Riesz space then the following lattice identities hold, *<sup>y</sup>* <sup>=</sup> *<sup>y</sup>*<sup>+</sup> <sup>−</sup> *<sup>y</sup>*<sup>−</sup> and *<sup>y</sup>*<sup>−</sup> = (−*y*)+. It is clear that *<sup>x</sup>* <sup>−</sup> *<sup>α</sup>***<sup>1</sup>** <sup>=</sup> (*<sup>x</sup>* <sup>−</sup> *<sup>α</sup>***1**)<sup>+</sup> <sup>−</sup> (*<sup>x</sup>* <sup>−</sup> *<sup>α</sup>***1**)<sup>−</sup> = (*<sup>x</sup>* <sup>−</sup> *<sup>α</sup>***1**)<sup>+</sup> <sup>−</sup> (*α***<sup>1</sup>** <sup>−</sup> *<sup>x</sup>*)<sup>+</sup> <sup>=</sup> *<sup>c</sup>*(*x*, *<sup>α</sup>*) <sup>−</sup> *<sup>p</sup>*(*x*, *<sup>α</sup>*). Therefore we have the identity

$$\mathbf{x} - a\mathbf{1} = c(\mathbf{x}, a) - p(\mathbf{x}, a)\mathbf{1}$$

which is called *put-call parity*.

If both *c*(*x*, *α*) > 0 and *p*(*x*, *α*) > 0, we say that the call option *c*(*x*, *α*) and the put option *p*(*x*, *α*) are *non trivial* and the exercise price *α* is a *non trivial exercise price* of *x*. If *c*(*x*, *α*) and *p*(*x*, *α*) belong to *X* then we say that *c*(*x*, *α*) and *p*(*x*, *α*) are *replicated*. If we suppose that **1** ∈ *X* and at least one of *c*(*x*, *a*), *p*(*x*, *a*) is replicated, then both of them are replicated since, *x* − *α***1** = *c*(*x*, *α*) − *p*(*x*, *α*).

#### **3.2. Completion of security markets**

We shall discuss the problem of completion by options of a two-period security market in which the space of marketed securities is a subspace of **R***m*. The present study involves vector sublattices generated by a subset *B* of **R***<sup>m</sup>* of positive, linearly independent vectors. A computational solution to this problem is provided by using the SUBlat Matlab function from [16].

#### 8 186 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 3 MATLAB Aided Option Replication <sup>9</sup>

Let us assume that in the beginning of a time period there are *n* securities traded in a market. Let <sup>S</sup> <sup>=</sup> {1, ..., *<sup>m</sup>*} denote a finite set of states and *xj* <sup>∈</sup> **<sup>R</sup>***<sup>m</sup>* <sup>+</sup> be the payoff vector of security *j* in *m* states. The payoffs *x*1, *x*2, ..., *xn* are assumed linearly independent so that there are no redundant securities. If *<sup>θ</sup>* = (*θ*1, *<sup>θ</sup>*2, ..., *<sup>θ</sup>n*) <sup>∈</sup> **<sup>R</sup>***<sup>n</sup>* is a non-zero portfolio then its payoff is the vector

**Example 2.** Suppose that in a security market, the payoff space is **R**<sup>12</sup> and the primitive

*x*<sup>1</sup> = (1, 2, 2, −1, 1, −2, −1, −3, 0, 0, 0, 0) *x*<sup>2</sup> = (0, 2, 0, 0, 1, 2, 0, 3, −1, −1, −1, −2) *x*<sup>3</sup> = (1, 2, 2, 0, 1, 0, 0, 0, −1, −1, −1, −2)

*u* = (1, 2, 2, 1, 1, 2, 1, 3, −1, −1, −1, −2).

where *X* denotes a matrix whose rows are the vectors *x*1, *x*2, *x*3, *u*. We can determine the completion by options of *X* i.e., the space *FU*(*X*), with the SUBlat function from [16] by

122010000000 020012030000 122112130000 000000001112 204000000000

000000001112 000100100000 000004060000 408000000000 060030000000

We consider a two-period security market *X* with a finite number *m* of states and a finite number of primitive securities with payoffs in **R***<sup>m</sup>* and we construct computational methods in order to determine maximal submarkets of *X* that replicate any option by using the results provided in subsection 2.2. In particular, in the theory of security markets it is a usual practice

**3.3. Computation of maximal submarkets that replicate any option**

<sup>1</sup> , *x*<sup>−</sup> <sup>1</sup> , *<sup>x</sup>*<sup>+</sup> <sup>2</sup> , *x*<sup>−</sup> <sup>2</sup> , *<sup>x</sup>*<sup>+</sup> <sup>3</sup> , *x*<sup>−</sup> <sup>3</sup> , *<sup>u</sup>*<sup>+</sup> <sup>1</sup> , *u*<sup>−</sup> 1 }

MATLAB Aided Option Replication 187

and that the strike subspace is the vector subspace *U* generated by the vector

>>XX = [max(X,zeros(size(X)));max(-X,zeros(size(X)))];

>>Linearindep = accumarray(I,J,[rank(XX),1],@min)';

>>[VectorSublattice,Positivebasis]=SUBlat(W')

Then, a maximal subset of linearly independent vectors of {*x*<sup>+</sup>

can be calculated by using the following code:

securities are:

>>S = rref(XX'); >>[I,J] = find(S);

using the following code:

VectorSublattice =

Positivebasis =

The results then are as follows

>>W = XX(Linearindep,:);

$$T(\theta) = \sum\_{i=1}^{n} \theta\_i \mathbf{x}\_i.$$

The set of payoffs of all portfolios is referred as the space of *marketed securities* and it is the linear span of the payoffs vectors *x*1, *x*2, ..., *xn* in **R***<sup>m</sup>* which we shall denote it by *X*, i.e.,

$$X = [\mathfrak{x}\_1, \mathfrak{x}\_2, \dots \mathfrak{x}\_n].$$

For any *<sup>x</sup>*, *<sup>u</sup>* <sup>∈</sup> **<sup>R</sup>***<sup>m</sup>* and any real number *<sup>a</sup>* the vector *cu*(*x*, *<sup>a</sup>*)=(*<sup>x</sup>* <sup>−</sup> *au*)<sup>+</sup> is the *call option* and *pu*(*x*, *<sup>a</sup>*)=(*au* <sup>−</sup> *<sup>x</sup>*)<sup>+</sup> is the *put option* of *<sup>x</sup>* with respect to the *strike vector u* and *exercise price a*.

Let *U* be a fixed subspace of **R***<sup>m</sup>* which is called *strike subspace* and the elements of *U* are the *strike vectors*. Then, the *completion by options* of the subspace *X* with respect to *U* is the space *FU*(*X*) which is defined inductively as follows:


$$\bullet \quad F\_{\mathcal{U}}(X) = \cup\_{n=1}^{\infty} X\_n.$$

The completion by options *FU*(*X*) of *X* with respect to *U* is the vector sublattice of **R***<sup>m</sup>* generated by the subspace *Y* = *X* ∪ *U*. The details are presented in the next theorem,

**Theorem 3.** *In the above notation, we have*

$$(i) \quad \mathcal{Y} \subseteq \mathcal{X}\_{1\prime}$$


Any set {*y*1, *<sup>y</sup>*2,..., *yr*} of linearly independent positive vectors of **<sup>R</sup>***<sup>m</sup>* such that *FU*(*X*) is the sublattice of **<sup>R</sup>***<sup>m</sup>* generated by {*y*1, *<sup>y</sup>*2,..., *yr*} is a *basic set* of the market.

**Theorem 4.** *Any maximal subset* {*y*1, *y*2,..., *yr*} *of linearly independent vectors of* A *is a basic set of the market, where* <sup>A</sup> <sup>=</sup> {*x*<sup>+</sup> <sup>1</sup> , *x*<sup>−</sup> <sup>1</sup> ,..., *<sup>x</sup>*<sup>+</sup> *<sup>n</sup>* , *x*<sup>−</sup> *<sup>n</sup>* }*, if U* ⊆ *X and* A = {*x*<sup>+</sup> <sup>1</sup> , *x*<sup>−</sup> <sup>1</sup> ,..., *<sup>x</sup>*<sup>+</sup> *<sup>n</sup>* , *x*<sup>−</sup> *<sup>n</sup>* , *<sup>u</sup>*<sup>+</sup> <sup>1</sup> , *u*<sup>−</sup> <sup>1</sup> ,..., *<sup>u</sup>*<sup>+</sup> *<sup>d</sup>* , *u*<sup>−</sup> *<sup>d</sup>* }*, if U X*

The space of marketed securities *X* is *complete by options* with respect to *U* if *X* = *FU*(*X*).

From theorem 1 it follows,

**Theorem 5.** *The space X of marketed securities is complete by options with respect to U if and only if U* ⊆ *X and cardR*(*β*) = *n. In addition, the dimension of FU*(*X*) *is equal to the cardinal number of R*(*β*)*. Therefore, FU*(*X*) = **R***<sup>k</sup> if and only if cardR*(*β*) = *k*.

**Example 2.** Suppose that in a security market, the payoff space is **R**<sup>12</sup> and the primitive securities are:

$$\begin{aligned} \mathbf{x}\_1 &= (1, 2, 2, -1, 1, -2, -1, -3, 0, 0, 0, 0) \\ \mathbf{x}\_2 &= (0, 2, 0, 0, 1, 2, 0, 3, -1, -1, -1, -2) \\ \mathbf{x}\_3 &= (1, 2, 2, 0, 1, 0, 0, 0, -1, -1, -1, -2) \end{aligned}$$

and that the strike subspace is the vector subspace *U* generated by the vector

$$
\mu = (1, 2, 2, 1, 1, 2, 1, 3, -1, -1, -1, -2).
$$

Then, a maximal subset of linearly independent vectors of {*x*<sup>+</sup> <sup>1</sup> , *x*<sup>−</sup> <sup>1</sup> , *<sup>x</sup>*<sup>+</sup> <sup>2</sup> , *x*<sup>−</sup> <sup>2</sup> , *<sup>x</sup>*<sup>+</sup> <sup>3</sup> , *x*<sup>−</sup> <sup>3</sup> , *<sup>u</sup>*<sup>+</sup> <sup>1</sup> , *u*<sup>−</sup> 1 } can be calculated by using the following code:

```
>>XX = [max(X,zeros(size(X)));max(-X,zeros(size(X)))];
>>S = rref(XX');
>>[I,J] = find(S);
>>Linearindep = accumarray(I,J,[rank(XX),1],@min)';
>>W = XX(Linearindep,:);
```
where *X* denotes a matrix whose rows are the vectors *x*1, *x*2, *x*3, *u*. We can determine the completion by options of *X* i.e., the space *FU*(*X*), with the SUBlat function from [16] by using the following code:

>>[VectorSublattice,Positivebasis]=SUBlat(W')

The results then are as follows

8

vector

• *FU*(*X*) = <sup>∪</sup><sup>∞</sup>

(*i*) *Y* ⊆ *X*1*,*

{*x*<sup>+</sup> <sup>1</sup> , *x*<sup>−</sup>

<sup>1</sup> ,..., *<sup>x</sup>*<sup>+</sup>

From theorem 1 it follows,

*<sup>n</sup>* , *x*<sup>−</sup> *<sup>n</sup>* , *<sup>u</sup>*<sup>+</sup> <sup>1</sup> , *u*<sup>−</sup>

Let us assume that in the beginning of a time period there are *n* securities traded in a market.

*j* in *m* states. The payoffs *x*1, *x*2, ..., *xn* are assumed linearly independent so that there are no redundant securities. If *<sup>θ</sup>* = (*θ*1, *<sup>θ</sup>*2, ..., *<sup>θ</sup>n*) <sup>∈</sup> **<sup>R</sup>***<sup>n</sup>* is a non-zero portfolio then its payoff is the

The set of payoffs of all portfolios is referred as the space of *marketed securities* and it is the linear span of the payoffs vectors *x*1, *x*2, ..., *xn* in **R***<sup>m</sup>* which we shall denote it by *X*, i.e.,

*X* = [*x*1, *x*2, ...*xn*].

For any *<sup>x</sup>*, *<sup>u</sup>* <sup>∈</sup> **<sup>R</sup>***<sup>m</sup>* and any real number *<sup>a</sup>* the vector *cu*(*x*, *<sup>a</sup>*)=(*<sup>x</sup>* <sup>−</sup> *au*)<sup>+</sup> is the *call option* and *pu*(*x*, *<sup>a</sup>*)=(*au* <sup>−</sup> *<sup>x</sup>*)<sup>+</sup> is the *put option* of *<sup>x</sup>* with respect to the *strike vector u* and *exercise price a*. Let *U* be a fixed subspace of **R***<sup>m</sup>* which is called *strike subspace* and the elements of *U* are the *strike vectors*. Then, the *completion by options* of the subspace *X* with respect to *U* is the space

• *<sup>X</sup>*<sup>1</sup> is the subspace of **<sup>R</sup>***<sup>m</sup>* generated by <sup>O</sup>1, where <sup>O</sup><sup>1</sup> <sup>=</sup> {*cu*(*x*, *<sup>a</sup>*)|*<sup>x</sup>* <sup>∈</sup> *<sup>X</sup>*, *<sup>u</sup>* <sup>∈</sup> *<sup>U</sup>*, *<sup>a</sup>* <sup>∈</sup> **<sup>R</sup>**},

• *Xn* is the subspace of **<sup>R</sup>***<sup>m</sup>* generated by <sup>O</sup>*n*, where <sup>O</sup>*<sup>n</sup>* <sup>=</sup> {*cu*(*x*, *<sup>a</sup>*)|*<sup>x</sup>* <sup>∈</sup> *Xn*−1, *<sup>u</sup>* <sup>∈</sup> *<sup>U</sup>*, *<sup>a</sup>* <sup>∈</sup> **<sup>R</sup>**},

The completion by options *FU*(*X*) of *X* with respect to *U* is the vector sublattice of **R***<sup>m</sup>*

Any set {*y*1, *<sup>y</sup>*2,..., *yr*} of linearly independent positive vectors of **<sup>R</sup>***<sup>m</sup>* such that *FU*(*X*) is the

**Theorem 4.** *Any maximal subset* {*y*1, *y*2,..., *yr*} *of linearly independent vectors of* A *is*

**Theorem 5.** *The space X of marketed securities is complete by options with respect to U if and only if U* ⊆ *X and cardR*(*β*) = *n. In addition, the dimension of FU*(*X*) *is equal to the cardinal number of*

*<sup>d</sup>* }*, if U X*

The space of marketed securities *X* is *complete by options* with respect to *U* if *X* = *FU*(*X*).

<sup>1</sup> , *x*<sup>−</sup>

<sup>1</sup> ,..., *<sup>x</sup>*<sup>+</sup>

*<sup>n</sup>* , *x*<sup>−</sup>

*<sup>n</sup>* }*, if U* ⊆ *X and* A =

generated by the subspace *Y* = *X* ∪ *U*. The details are presented in the next theorem,

*n* ∑ *i*=1 *θixi*.

*T*(*θ*) =

<sup>+</sup> be the payoff vector of security

Let <sup>S</sup> <sup>=</sup> {1, ..., *<sup>m</sup>*} denote a finite set of states and *xj* <sup>∈</sup> **<sup>R</sup>***<sup>m</sup>*

*FU*(*X*) which is defined inductively as follows:

*<sup>n</sup>*=1*Xn*.

**Theorem 3.** *In the above notation, we have*

(*ii*) *FU*(*X*) *is the sublattice S*(*Y*) *of* **R***<sup>m</sup> generated by Y, and*

*a basic set of the market, where* <sup>A</sup> <sup>=</sup> {*x*<sup>+</sup>

(*iii*) *if U* <sup>⊆</sup> *X, then FU*(*X*) *is the sublattice of* **<sup>R</sup>***<sup>m</sup> generated by X*.

<sup>1</sup> ,..., *<sup>u</sup>*<sup>+</sup>

*R*(*β*)*. Therefore, FU*(*X*) = **R***<sup>k</sup> if and only if cardR*(*β*) = *k*.

sublattice of **<sup>R</sup>***<sup>m</sup>* generated by {*y*1, *<sup>y</sup>*2,..., *yr*} is a *basic set* of the market.

*<sup>d</sup>* , *u*<sup>−</sup>

denotes the set of call options written on the elements of *X*,

denotes the set of call options written on the elements of *Xn*−1,


#### **3.3. Computation of maximal submarkets that replicate any option**

We consider a two-period security market *X* with a finite number *m* of states and a finite number of primitive securities with payoffs in **R***<sup>m</sup>* and we construct computational methods in order to determine maximal submarkets of *X* that replicate any option by using the results provided in subsection 2.2. In particular, in the theory of security markets it is a usual practice to take call and put options with respect to the riskless bond **1** = (1, 1, ..., 1). Then, the completion *F***1**(*X*) of *X* by options (see subsection 3.2) is the subspace of **R***<sup>m</sup>* generated by all options written on the elements of *<sup>X</sup>* ∪ {**1**}. Since the payoff space is **<sup>R</sup>***m*, which is a vector lattice, in the case where **1** ∈ *X* then *F***1**(*X*) is exactly the vector sublattice generated by *X*. If, in addition, *X* is a vector sublattice of **R***<sup>m</sup>* then *F***1**(*X*) = *X* therefore any option is replicated.

**Theorem 6.** *Let* {*bi*, *i* = 1, 2, ..., *μ*} *be the positive basis of F***1**(*X*) *which is a partition of the unit and*

*(ii) there exists a maximal proper partition δ* = {*σi*|*i* = 1, 2, ..., *k*} *of* {1, 2, ..., *n*} *so that Y is the*

We present the proposed computational method that enables us to determine maximal submarkets that replicate any option. Our numerical method is based on the introduction of the mrsubspace function, from [15], that allow us to perform fast testing for a variety of

Recall that *X* is the security market generated by a collection {*x*1, *x*2, ..., *xn*} of linearly independent vectors (not necessarily positive) of **<sup>R</sup>***m*. If **<sup>1</sup>** <sup>∈</sup> *<sup>X</sup>* then it is possible to determine a basic set of marketed securities i.e., a set of linearly independent and positive vectors of *X*.

**Proposition 3.** *If a* = max{�*xi*�∞|*i* = 1, 2, ..., *n*}*, then at least one of the two sets of positive vectors*

{*yi* = *a***1** − *xi*|*i* = 1, 2, ..., *n*}, {*zi* = 2*a***1** − *xi*|*i* = 1, 2, ..., *n*},

The main steps of the underlying algorithmic procedure that enables us to determine maximal

(1) Use proposition 3 in order to determine a basic set {*y*1, *y*2, ..., *yn*} of marketed securities.

(4) Use Theorem 1 in order to construct the vector sublattice generated by *y*1, *y*2, ..., *yn*, which is exactly the completion by options *F***1**(*X*) of *X*. Then, determine a positive basis

(6) Use Proposition 1 in order to determine a positive basis {*b*1, *b*2, ..., *bμ*} of *F***1**(*X*) which is a

(8) Decide which of the proper partitions created in step (7) are maximal proper partitions

 1, *d* 2, ..., *d*

> 1, *b* 2, ..., *b n*}.

*<sup>n</sup>*} of *X*.

MATLAB Aided Option Replication 189

(2) Use Equation (1) in order to determine the basic curve *β* of the vectors *yi*.

**4.1. Algorithm for calculating maximal submarkets that replicate any option**

*<sup>i</sup>*, *i* = 1, 2, ..., *n*} *be the projection basis of X corresponding to the basis* {*bi*}*. If Y is a subspace of*

*let* {*b* 

*of X*

*X, the following are equivalent:*

*(i) Y is a maximal replicated subspace of X,*

*The set of maximal replicated submarkets of X is nonempty.*

This is possible through the following easy proposition:

submarkets that replicate any option are the following:

(5) Use Theorem 2 in order to determine a projection basis {*d*

partition of the unit and the corresponding projection basis {*b*

(7) Calculate all the possible proper partitions of the set {1, 2, ..., *n*}.

and determine the corresponding maximal replicated submarkets.

*sublattice of* **R***<sup>m</sup> generated by δ.*

**4. The computational method**

*consists of linearly independent vectors.*

(3) Determine the range *R*(*β*) of *β*.

{*d*1, *d*2, ..., *dμ*} for *F***1**(*X*).

dimensions and subspaces.

A basic set of marketed securities (i.e., a set of linearly independent and positive vectors) of *X* always exist and the sublattice of **R***<sup>m</sup>* generated by a basic set of marketed securities is *F***1**(*X*). In addition, *F***1**(*X*) has a positive basis which is a partition of the unit.

Let us assume that *X* is generated by a basic set of marketed securities, then from Theorem 1 it is possible to determine a positive basis {*b*1, *b*2, ..., *bμ*} of *F***1**(*X*).

The sublattice *Z*, generated by a basic set of marketed securities, is exactly *F***1**(*X*) and *F***1**(*X*) has a positive basis which is a partition of the unit, i.e., <sup>∑</sup>*<sup>μ</sup> <sup>i</sup>*=<sup>1</sup> *bi* = **1**. This is possible since the notion of a positive basis is unique in the sense of positive multiples therefore we are able to extract from the positive basis {*b*1, *b*2, ..., *bμ*} another positive basis {*d*1, *d*2, ..., *dμ*} of *F***1**(*X*) which is a partition of the unit. Therefore, let us denote by {*d*1, *d*2, ..., *dμ*} a positive basis of *F***1**(*X*) which is a partition of the unit. Then, by Theorem 2, if

$$(\tilde{d}\_1, \tilde{d}\_2, \dots, \tilde{d}\_r)^T = A^{-1} (z\_1, z\_2, \dots, z\_r)^T \dots$$

where *A* is the *r* × *r* matrix whose *i*th column is the vector *Pi*, for *i* = 1, 2, ...,*r* then {*d* 1, *d* 2, ..., *d <sup>r</sup>*} is a projection basis of *F***1**(*X*). The projection basis {*d* 1, *d* 2, ..., *d <sup>r</sup>*} is the projection basis of *X* corresponding to the basis {*d*1, *d*2, ..., *dμ*}. For *Z* = *F***1**(*X*) proposition 1 takes the following form.

**Proposition 2.** *Suppose that* {*di*} *is the basis of F***1**(*X*) *given by equation (2) of theorem 1 and* {*d i*} *is the projection basis of X corresponding to the basis* {*di*}. *Then* {*bi* <sup>=</sup> *di* �*di*�<sup>∞</sup> <sup>|</sup>*<sup>i</sup>* <sup>=</sup> 1, 2, ..., *<sup>μ</sup>*} *is the positive basis of F***1**(*X*) *which is a partition of the unit and* {*b <sup>i</sup>* = *<sup>d</sup> i* �*di*�<sup>∞</sup> <sup>|</sup>*<sup>i</sup>* <sup>=</sup> 1, 2, ..., *<sup>n</sup>*} *is the projection basis of X corresponding to the basis* {*bi*} *of F***1**(*X*)*.*

Suppose that *Y* is a subspace of *X*, then if *F***1**(*Y*) ⊆ *X* we say that *Y* is *replicated*. If, in addition, for any subspace *Z* of *X* with *Y Z* we have that *X F***1**(*Z*) then *Y* is a *maximal replicated subspace* or a *maximal replicated submarket* of *X*. Note that, the *replicated kernel of the market*, i.e., the union of all maximal replicated subspaces of the market is the set of all elements *x* of *X* so that any option written on *x* is replicated.

Let {*bi i* = 1, 2, ..., *μ*} be a positive basis of *F***1**(*X*) which is a partition of the unit and let {*b <sup>i</sup> i* = 1, 2, ..., *n*} be the projection basis of *X* corresponding to the basis {*bi*}. Recall that, a partition *δ* = {*σi*|*i* = 1, 2, ..., *k*} of {1, 2, ..., *n*} is *proper* if for any *r* = 1, 2, ..., *k*, the vector *wr* = <sup>∑</sup>*i*∈*σ<sup>r</sup> <sup>b</sup> <sup>i</sup>* is a binary vector with ∑*<sup>k</sup> <sup>r</sup>*=<sup>1</sup> *wr* = **1**. If there is no proper partition of {1, 2, ..., *n*} strictly finer than *δ*, then we say that *δ* is a *maximal proper partition* of {1, 2, ..., *n*}.

The following theorem provides the development of a method in order to determine the set of securities with replicated options by using the theory of positive bases and projection bases.

**Theorem 6.** *Let* {*bi*, *i* = 1, 2, ..., *μ*} *be the positive basis of F***1**(*X*) *which is a partition of the unit and let* {*b <sup>i</sup>*, *i* = 1, 2, ..., *n*} *be the projection basis of X corresponding to the basis* {*bi*}*. If Y is a subspace of X, the following are equivalent:*

*(i) Y is a maximal replicated subspace of X,*

10

{*d* 1, *d* 2, ..., *d*

{*b* 

*wr* = <sup>∑</sup>*i*∈*σ<sup>r</sup> <sup>b</sup>*

following form.

to take call and put options with respect to the riskless bond **1** = (1, 1, ..., 1). Then, the completion *F***1**(*X*) of *X* by options (see subsection 3.2) is the subspace of **R***<sup>m</sup>* generated by all options written on the elements of *<sup>X</sup>* ∪ {**1**}. Since the payoff space is **<sup>R</sup>***m*, which is a vector lattice, in the case where **1** ∈ *X* then *F***1**(*X*) is exactly the vector sublattice generated by *X*. If, in addition, *X* is a vector sublattice of **R***<sup>m</sup>* then *F***1**(*X*) = *X* therefore any option is replicated. A basic set of marketed securities (i.e., a set of linearly independent and positive vectors) of *X* always exist and the sublattice of **R***<sup>m</sup>* generated by a basic set of marketed securities is *F***1**(*X*).

Let us assume that *X* is generated by a basic set of marketed securities, then from Theorem 1

The sublattice *Z*, generated by a basic set of marketed securities, is exactly *F***1**(*X*) and *F***1**(*X*)

the notion of a positive basis is unique in the sense of positive multiples therefore we are able to extract from the positive basis {*b*1, *b*2, ..., *bμ*} another positive basis {*d*1, *d*2, ..., *dμ*} of *F***1**(*X*) which is a partition of the unit. Therefore, let us denote by {*d*1, *d*2, ..., *dμ*} a positive basis of

where *A* is the *r* × *r* matrix whose *i*th column is the vector *Pi*, for *i* = 1, 2, ...,*r* then

basis of *X* corresponding to the basis {*d*1, *d*2, ..., *dμ*}. For *Z* = *F***1**(*X*) proposition 1 takes the

**Proposition 2.** *Suppose that* {*di*} *is the basis of F***1**(*X*) *given by equation (2) of theorem 1 and* {*d*

Suppose that *Y* is a subspace of *X*, then if *F***1**(*Y*) ⊆ *X* we say that *Y* is *replicated*. If, in addition, for any subspace *Z* of *X* with *Y Z* we have that *X F***1**(*Z*) then *Y* is a *maximal replicated subspace* or a *maximal replicated submarket* of *X*. Note that, the *replicated kernel of the market*, i.e., the union of all maximal replicated subspaces of the market is the set of all elements *x* of *X* so

Let {*bi i* = 1, 2, ..., *μ*} be a positive basis of *F***1**(*X*) which is a partition of the unit and let

The following theorem provides the development of a method in order to determine the set of securities with replicated options by using the theory of positive bases and projection bases.

strictly finer than *δ*, then we say that *δ* is a *maximal proper partition* of {1, 2, ..., *n*}.

*<sup>i</sup> i* = 1, 2, ..., *n*} be the projection basis of *X* corresponding to the basis {*bi*}. Recall that, a partition *δ* = {*σi*|*i* = 1, 2, ..., *k*} of {1, 2, ..., *n*} is *proper* if for any *r* = 1, 2, ..., *k*, the vector

*<sup>r</sup>*} is a projection basis of *F***1**(*X*). The projection basis {*d*

*is the projection basis of X corresponding to the basis* {*di*}. *Then* {*bi* <sup>=</sup> *di*

*<sup>r</sup>*)*<sup>T</sup>* <sup>=</sup> *<sup>A</sup>*−1(*z*1, *<sup>z</sup>*2, ..., *zr*)*T*,

 *<sup>i</sup>* = *<sup>d</sup> i*

*<sup>i</sup>*=<sup>1</sup> *bi* = **1**. This is possible since

*<sup>r</sup>*} is the projection

�*di*�<sup>∞</sup> <sup>|</sup>*<sup>i</sup>* <sup>=</sup> 1, 2, ..., *<sup>μ</sup>*} *is the*

�*di*�<sup>∞</sup> <sup>|</sup>*<sup>i</sup>* <sup>=</sup> 1, 2, ..., *<sup>n</sup>*} *is the projection*

 *i*}

 1, *d* 2, ..., *d*

*<sup>r</sup>*=<sup>1</sup> *wr* = **1**. If there is no proper partition of {1, 2, ..., *n*}

In addition, *F***1**(*X*) has a positive basis which is a partition of the unit.

it is possible to determine a positive basis {*b*1, *b*2, ..., *bμ*} of *F***1**(*X*).

has a positive basis which is a partition of the unit, i.e., <sup>∑</sup>*<sup>μ</sup>*

*F***1**(*X*) which is a partition of the unit. Then, by Theorem 2, if

(*d* 1, *d* 2, ..., *d*

*positive basis of F***1**(*X*) *which is a partition of the unit and* {*b*

*basis of X corresponding to the basis* {*bi*} *of F***1**(*X*)*.*

that any option written on *x* is replicated.

*<sup>i</sup>* is a binary vector with ∑*<sup>k</sup>*

*(ii) there exists a maximal proper partition δ* = {*σi*|*i* = 1, 2, ..., *k*} *of* {1, 2, ..., *n*} *so that Y is the sublattice of* **R***<sup>m</sup> generated by δ.*

*The set of maximal replicated submarkets of X is nonempty.*

#### **4. The computational method**

We present the proposed computational method that enables us to determine maximal submarkets that replicate any option. Our numerical method is based on the introduction of the mrsubspace function, from [15], that allow us to perform fast testing for a variety of dimensions and subspaces.

#### **4.1. Algorithm for calculating maximal submarkets that replicate any option**

Recall that *X* is the security market generated by a collection {*x*1, *x*2, ..., *xn*} of linearly independent vectors (not necessarily positive) of **<sup>R</sup>***m*. If **<sup>1</sup>** <sup>∈</sup> *<sup>X</sup>* then it is possible to determine a basic set of marketed securities i.e., a set of linearly independent and positive vectors of *X*. This is possible through the following easy proposition:

**Proposition 3.** *If a* = max{�*xi*�∞|*i* = 1, 2, ..., *n*}*, then at least one of the two sets of positive vectors of X*

$$\{y\_i = a\mathbf{1} - \mathbf{x}\_i | i = 1, 2, \dots, n\}, \; \{z\_i = 2a\mathbf{1} - \mathbf{x}\_i | i = 1, 2, \dots, n\},$$

*consists of linearly independent vectors.*

The main steps of the underlying algorithmic procedure that enables us to determine maximal submarkets that replicate any option are the following:


In [15], it is presented the translation followed by the implementation of the aforementioned algorithm in **R***<sup>m</sup>* within a Matlab-based function named mrsubspace. Moreover, in the same paper, computational experiments assessing the effectiveness of this function and lead us to the conclusion that the mrsubspace function provides an important tool in order to investigate replicated subspaces.

**Example 3.** Consider the following 5 vectors *x*1, *x*2, ..., *x*<sup>5</sup> in **R**10,

⎡ ⎢ ⎢ ⎢ ⎢ ⎣

above collection of vectors we use the following simple code:

as a result and after removing irrelevant Matlab output one gets

>> X = [0,1,0,1,1,1,1,2,2,1;1,1,1,2,1,1,1,2,1,2;...

1,1,1,2,1,1,1,2,1,1;1,1,1,1,1,1,1,2,2,1;2,1,2,1,1,1,1,1,1,1];

1010000000 0101111000 0000000110 0000000001

1010000000 0100111010 0001000100 0000000001

1010000000 0100111000 0001000000 0000000100 0000000001 0000000010

and *X* = [*x*1, *x*2, ..., *x*5] is the marketed space.

>> [Npb,Cprb] =mrsubspace(X)

The 1 partition(s) are: {1} {2 3} {4} {5}

The 1 partition(s) are: {1} {2} {3 4} {5}

ReplicatedSubspace =

Npb =

ReplicatedSubspace =

*x*1 *x*2 *x*3 *x*4 *x*5 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣

0101111221 1112111212 1112111211 1111111221 2121111111

Note that **1** = *x*<sup>5</sup> − *x*<sup>4</sup> + *x*1. In order to determine the maximal replicated subspaces for the

⎤ ⎥ ⎥ ⎥ ⎥ ⎦ MATLAB Aided Option Replication 191

## **4.2. The computational approach - Code features**

We shall discuss the proposed computational method that enables us to determine maximal submarkets that replicate any option. The standard method used currently to determine the maximal replicated submarkets is based on a manual processing. From section 3, it is evident that the required number of verifications for this process can be of significant size even in a relatively low-dimensional space, thus rendering the problem too difficult to solve. Our numerical method is based on the introduction of the mrsubspace function, from [15] that allow us to perform fast testing for a variety of dimensions and subspaces.

The structure of the code ensures flexibility, meaning that it is convenient for applications as well as for research and educational purposes. The given security market *X*, generated by the linearly independent vectors *x*1, *x*2, ..., *xn*, must be given under a matrix notation with columns the vectors *x*1, *x*2, ..., *xn*. The mrsubspace function must be stored in a Matlab-accessible directory and then the input data, i.e., the matrix *X*, can be typed directly in the Matlab's environment. Under the following command,

mrsubspace(X);

the program solves the problem of option replication and prints out the maximal proper partitions as well as the corresponding maximal replicated subspaces. If *X* is a vector sublattice, then *X* = *F***1**(*X*) and any option is replicated. In the case where the initial space *X* is not a vector sublattice, it is possible to produce the normalized positive basis and the corresponding projection basis with the following code,

```
[Npb,Cprb] = mrsubspace(X)
```
Inside the code there are several explanations that indicate the implemented part of the algorithm. A user proficient in Matlab can easily use the code and modify it if needed. Especially, the user can isolate a part of the code according to his/her special needs to solve different problems like


In the last part of the code, entitled Maximal proper partitions - Maximal replicated subspaces, the user can change the way that the mrsubspace function understands the values 0 and 1, according to his/her knowledge and needs.

**Example 3.** Consider the following 5 vectors *x*1, *x*2, ..., *x*<sup>5</sup> in **R**10,

$$
\begin{bmatrix}
\mathbf{x}\_1 \\
\mathbf{x}\_2 \\
\mathbf{x}\_3 \\
\mathbf{x}\_4 \\
\mathbf{x}\_5
\end{bmatrix} = \begin{bmatrix}
0 \ 1 \ 0 \ 1 \ 1 \ 1 \ 1 \ 2 \ 2 \ 1 \\
1 \ 1 \ 1 \ 2 \ 1 \ 1 \ 1 \ 2 \ 1 \ 2 \\
1 \ 1 \ 1 \ 2 \ 1 \ 1 \ 1 \ 2 \ 1 \ 1 \\
1 \ 1 \ 1 \ 1 \ 1 \ 1 \ 2 \ 2 \ 1 \\
2 \ 1 \ 2 \ 1 \ 1 \ 1 \ 1 \ 1 \ 1 \ 1 \ 1
\end{bmatrix}
$$

and *X* = [*x*1, *x*2, ..., *x*5] is the marketed space.

12

investigate replicated subspaces.

**4.2. The computational approach - Code features**

environment. Under the following command,

[Npb,Cprb] = mrsubspace(X)

different problems like

corresponding projection basis with the following code,

• Determine a basic set of marketed securities.

according to his/her knowledge and needs.

by a finite collection of linearly independent vectors of **R***m*.

mrsubspace(X);

In [15], it is presented the translation followed by the implementation of the aforementioned algorithm in **R***<sup>m</sup>* within a Matlab-based function named mrsubspace. Moreover, in the same paper, computational experiments assessing the effectiveness of this function and lead us to the conclusion that the mrsubspace function provides an important tool in order to

We shall discuss the proposed computational method that enables us to determine maximal submarkets that replicate any option. The standard method used currently to determine the maximal replicated submarkets is based on a manual processing. From section 3, it is evident that the required number of verifications for this process can be of significant size even in a relatively low-dimensional space, thus rendering the problem too difficult to solve. Our numerical method is based on the introduction of the mrsubspace function, from [15] that

The structure of the code ensures flexibility, meaning that it is convenient for applications as well as for research and educational purposes. The given security market *X*, generated by the linearly independent vectors *x*1, *x*2, ..., *xn*, must be given under a matrix notation with columns the vectors *x*1, *x*2, ..., *xn*. The mrsubspace function must be stored in a Matlab-accessible directory and then the input data, i.e., the matrix *X*, can be typed directly in the Matlab's

the program solves the problem of option replication and prints out the maximal proper partitions as well as the corresponding maximal replicated subspaces. If *X* is a vector sublattice, then *X* = *F***1**(*X*) and any option is replicated. In the case where the initial space *X* is not a vector sublattice, it is possible to produce the normalized positive basis and the

Inside the code there are several explanations that indicate the implemented part of the algorithm. A user proficient in Matlab can easily use the code and modify it if needed. Especially, the user can isolate a part of the code according to his/her special needs to solve

• Find the completion *F***1**(*X*) of *X* by options in **R***<sup>m</sup>* or find the vector sublattice generated

• Calculate a positive basis and a projection basis for a finite dimensional vector sublattice.

In the last part of the code, entitled Maximal proper partitions - Maximal replicated subspaces, the user can change the way that the mrsubspace function understands the values 0 and 1,

allow us to perform fast testing for a variety of dimensions and subspaces.

Note that **1** = *x*<sup>5</sup> − *x*<sup>4</sup> + *x*1. In order to determine the maximal replicated subspaces for the above collection of vectors we use the following simple code:

>> X = [0,1,0,1,1,1,1,2,2,1;1,1,1,2,1,1,1,2,1,2;... 1,1,1,2,1,1,1,2,1,1;1,1,1,1,1,1,1,2,2,1;2,1,2,1,1,1,1,1,1,1]; >> [Npb,Cprb] =mrsubspace(X)

as a result and after removing irrelevant Matlab output one gets

```
The 1 partition(s) are:
 {1} {2 3} {4} {5}
ReplicatedSubspace =
 1010000000
 0101111000
 0000000110
 0000000001
The 1 partition(s) are:
 {1} {2} {3 4} {5}
ReplicatedSubspace =
 1010000000
 0100111010
 0001000100
 0000000001
Npb =
 1010000000
 0100111000
 0001000000
 0000000100
 0000000001
 0000000010
```
14 192 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 3 MATLAB Aided Option Replication <sup>15</sup>


Therefore, the marketed space *X* has two maximal replicated subspaces, {1} {2 3} {4} {5} and {1} {2} {3 4} {5} are maximal proper partitions with corresponding maximal replicated subspaces the subspaces

**6. References**

[1] Aliprantis, C.D. & Brown, D.J. (1983). Equilibria in markets with a Riesz space of

MATLAB Aided Option Replication 193

[2] Aliprantis, C.D., Brown, D.J. & Burkinshaw, O. (1987a). Edgeworth equilibria,

[3] Aliprantis, C.D., Brown, D.J. & Burkinshaw, O. (1987b). Edgeworth equilibria in

[4] Aliprantis, C.D., Brown, D.J. & Burkinshaw, O. (1990). Existence and optimally of

[5] Aliprantis, C.D. & Burkinshaw, O. (1991). When is the core equivalence theorem valid?,

[6] Aliprantis, C.D., Tourky, R. & Yannelis, N. C. (2001). A theory of value with non–linear prices. Equilibrium analysis beyond vector lattices, *J. Econom. Theory*, 100, pp. 22–72. [7] Aliprantis, C.D. & Tourky, R. (2002). Markets that don't replicate any option, *Economics*

[8] Baptista A.M. (2003). Spanning with american options, *Journal of Economic Theory*, 110,

[9] Baptista A.M. (2005). Options and efficiency in multidate security markets, *Mathematical*

[10] Baptista A.M. (2007). On the non-existence of redundant options, *Economic Theory*, 31,

[11] Katsikis, V.N. (2007). Computational methods in portfolio insurance, *Applied*

[12] Katsikis, V.N. (2008). Computational methods in lattice-subspaces of *C*[*a*, *b*] with applications in portfolio insurance, *Applied Mathematics and Computation*, 200,

[13] Katsikis, V.N. (2009). A Matlab-based rapid method for computing lattice-subspaces and vector sublattices of **R***n*: Applications in portfolio insurance, *Applied Mathematics*

[14] Katsikis, V.N. (2010). Computational and Mathematical Methods in Portfolio Insurance. A MATLAB-Based Approach., Matlab - Modelling, Programming and Simulations, Emilson Pereira Leite(Ed.), ISBN: 978-953-307-125-1, InTech, Available from: http://www.intechopen.com/articles/show/title/computational-and-

[15] Katsikis, V.N. (2011). Computational methods for option replication, *International*

[16] Katsikis V.N. & Polyrakis I. (2012). Computation of vector sublattices and minimal lattice-subspaces. Applications in finance, *Applied Mathematics and Computation*, 218,

[17] Mas-Colell, A. (1986). The price equilibrium existence problem in topological vector

[18] Mas-Colell, A. & Richard, S.F. (1991). A new approach to the existence of equilibrium in

[19] Podczeck, K. (1996). Equilibria in vector lattices without ordered preferences or uniform

mathematical-methods-in-portfolio-insurance-a-matlab-based-approach-

*Journal of Computer Mathematics*, 88, No. 13, pp. 2752-2769.

commodities, *J. Math. Econom.*, 11, pp. 189–207.

production economies, *J. Econom. Theory*, 43, pp. 253–291.

*Econometrica*, 55, pp. 1109–1137.

*Econom. Theory*, 1, pp. 169–182.

*Fi- nance*, 15, No.4, pp.569-587.

*and Computation*, 215, pp.961-972.

lattices, *Econometrica*, 54, pp. 1039–1053.

vector lattices, *J. Econom. Theory*, 53, pp. 1–11.

properness, *J. Math. Econom.*, 25, pp. 465–484.

*Mathematics and Computation*, 189, pp.9-22.

*Letter*, 76, pp.443–447.

pp.264-289.

pp.205-212.

pp.204-219.

pp.6860–6873.

competitive equilibria, *Springer–Verlag*.

$$Y\_1 = [(1,0,1,0,0,0,0,0,0,0), (0,1,0,1,1,1,1,0,0,0), (0,0,0,0,0,0,0,1,1,0), \\
\begin{aligned} &(0,0,0,0,0,0,0,0,1) \\ &(1,0,0,0,0,0,1) \end{aligned}$$

and

$$Y\_2 = [(1,0,1,0,0,0,0,0,0,0), (0,1,0,0,1,1,1,0,1,0), (0,0,0,1,0,0,0,1,0,0), \\
\begin{aligned} &(0,0,0,0,0,0,0,0,1), \\ &(0,0,0,0,0,0,0,1), \end{aligned} \end{aligned}$$

respectively. The replicated kernel of the market is *Y* = *Y*<sup>1</sup> ∪ *Y*2.

## **5. Conclusions**

In this chapter, computational methods for option replication are presented based on vector lattice theory. It is well accepted that the lattice theoretic ideas are one of the most important technical contributions of the large literature on modern mathematical finance. In this chapter, we consider an incomplete market of primitive securities, meaning that some call and put options need not be marketed. Our objective is to describe an efficient method for computing maximal submarkets that replicate any option. Even though, there are several important results on option replication they cannot provide a method for the determination of the replicated options. By using the theory of vector-lattices and positive bases it is provided a procedure in order to determine the set of securities with replicated options. Moreover, we determine those subspaces of the marketed subspace that replicate any option by introducing a Matlab function, namely mrsubspace. The results of this work can give us an important tool in order to study the interesting problem of option replication of a two-period security market in which the space of marketed securities is a subspace of **R***m*. This work is based on a recent work, [15], regarding computational methods for option replication.

## **Author details**

Vasilios N. Katsikis *General Department of Mathematics, Technological Education Institute of Piraeus, 12244 Athens, Greece*

#### **6. References**

14

Cprb =

and

subspaces the subspaces

**5. Conclusions**

**Author details** Vasilios N. Katsikis

*Greece*

1010000000 0100111010 0 0 0 1 0 0 0 0 -1 0 0000000110 0000000001

Therefore, the marketed space *X* has two maximal replicated subspaces, {1} {2 3} {4} {5} and {1} {2} {3 4} {5} are maximal proper partitions with corresponding maximal replicated

> *Y*<sup>1</sup> = [(1, 0, 1, 0, 0, 0, 0, 0, 0, 0),(0, 1, 0, 1, 1, 1, 1, 0, 0, 0),(0, 0, 0, 0, 0, 0, 0, 1, 1, 0), (0, 0, 0, 0, 0, 0, 0, 0, 0, 1)]

> *Y*<sup>2</sup> = [(1, 0, 1, 0, 0, 0, 0, 0, 0, 0),(0, 1, 0, 0, 1, 1, 1, 0, 1, 0),(0, 0, 0, 1, 0, 0, 0, 1, 0, 0), (0, 0, 0, 0, 0, 0, 0, 0, 0, 1)],

In this chapter, computational methods for option replication are presented based on vector lattice theory. It is well accepted that the lattice theoretic ideas are one of the most important technical contributions of the large literature on modern mathematical finance. In this chapter, we consider an incomplete market of primitive securities, meaning that some call and put options need not be marketed. Our objective is to describe an efficient method for computing maximal submarkets that replicate any option. Even though, there are several important results on option replication they cannot provide a method for the determination of the replicated options. By using the theory of vector-lattices and positive bases it is provided a procedure in order to determine the set of securities with replicated options. Moreover, we determine those subspaces of the marketed subspace that replicate any option by introducing a Matlab function, namely mrsubspace. The results of this work can give us an important tool in order to study the interesting problem of option replication of a two-period security market in which the space of marketed securities is a subspace of **R***m*. This work is based on

*General Department of Mathematics, Technological Education Institute of Piraeus, 12244 Athens,*

a recent work, [15], regarding computational methods for option replication.

respectively. The replicated kernel of the market is *Y* = *Y*<sup>1</sup> ∪ *Y*2.

	- [20] Richard, S.F. (1989). A new approach to production equilibria in vector lattices, *J. Math. Econom.*, 18, pp. 231–247.
	- [21] Ross, S.A. (1976). Options and efficiency, *Quaterly Journal of Economics*, 90, pp.75-89.
	- [22] Tourky, R. (1998). A new approach to the limit theorem on the core of an economy in vector lattices, *J. Econom. Theory*, 78, pp. 321–328.

**Convolution Kernel for Fast CPU/GPU**

Sébastien Leclaire, Maud El-Hachem and Marcelo Reggio

Additional information is available at the end of the chapter

definition in the modeling of multiphase flows.

cited.

http://dx.doi.org/10.5772/46455

**1. Introduction**

**on a Square/Cubic Lattice**

**Computation of 2D/3D Isotropic Gradients**

**Chapter 9**

The design of discrete operators or filters for the calculation of gradients is a classical topic in scientific computing. Typical applications are gradient reconstruction in computational fluid dynamics, edge detection in computer graphics and biomedical imaging, and phase boundary

Edge detection, which is widely performed in image analysis, is an operation that requires gradient calculation. Commonly used edge detection methods are *Canny, Prewitt, Roberts* and *Sobel*, which can be found in MATLAB's platform. In this field, edge detection techniques rely on the application of convolution masks to provide a filter or kernel to calculate gradients in

For multiphase flows, an edge or contour corresponds to the interface between the fluids. In this respect, traditional gradient calculation methods based on 1D edge detection are not necessarily suited for the underlying physics, because there is no direction in which the gradients of the phase contours tend to evolve over time. As a result, definition of the geometric progress of the interface requires many gradient estimation computations, as is the case in moving and deforming bubbles or droplets, for example. Although it can still be a viable tool, it is clear that the 1D-based method is becoming less useful for simulating these

To address this issue, we present an efficient computational method for obtaining discrete isotropic gradients that was previously applied to simulate two-phase flows within the lattice Boltzman framework [1, 2]. This "omnidirectional" approach makes it possible to improve the limitations inherent in handling high density ratios between the phases and to significantly reduce spurious currents at the interface. The method is based on a filter which is generally not split along any direction, and there is no need to make the assumption of a continuous filter to reach isotropy, as done by [3]. We also believe that optimal or maximal isotropy can

> ©2012 Reggio et al., licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0),which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly

©2012 Reggio et al., licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

two perpendicular directions. A threshold is then applied to obtain an edge shape.

phenomena, which are not, in general, biased toward any particular direction.

[23] Yannelis, N.C. & Zame, W.R. (1986). Equilibria in Banach lattices without ordered preferences, *J. Math. Econom.*, 15, pp. 85–110.

194 MATLAB – A Fundamental Tool for Scienti c Computing and Engineering Applications – Volume 3 **Chapter 0 Chapter 9**
