Ane models in financial markets : pricing options using inversion techniques and the properties of affine models

(1)

Affine models in financial markets:

Pricing options using inversion techniques

and the properties of affine models

Frank Pardoel

July 24, 2016

Master Thesis Actuarial Science and Mathematical Finance Supervisor: prof. dr. F. Kleibergen

Amsterdam Business School Amsterdam Executive Programme

(2)

Statement of Originality

This document is written by Frank Pardoel who declares to take full responsibility for the contents of this document. I declare that the text and the work presented in this document is original and that no sources other than those mentioned in the text and its references have been used in creating it. The Faculty of Economics and Business is responsible solely for the supervision of completion of the work, not for the contents.

(3)

1. Introduction

An option gives the holder the right to buy or sell an underlying asset for a certain price at a certain moment in time. In financial markets, a variety of options are traded, e.g. options on stocks, options on futures, options on swaps. In order to determine the theoretical, unique price of a certain option, the no-arbitrage condition is used. The no-arbitrage condition formulates that one cannot earn a non-negative payoff in the fu-ture without investing a certain amount today. The fundamental no-arbitrage concept in relation to option pricing has been introduced in financial markets by Black and Sc-holes (1973). Together with the work of Merton (1973), this resulted in the Nobel Prize winning Black, Scholes and Merton option pricing model.

The great advantage of the Black, Scholes and Merton option pricing model is the analytical solution. Given the model parameters, one can directly determine the price of the option. The cost of the analytical traceability are the simplified underlying assump-tions, e.g. lognormal stock prices, constant volatility for each time and strike, constant interest rates. To relax the underlying assumptions, different alternative approaches have been explored. One alternative approach applied in financial mathematics is the use of binomial trees. Binomial trees use the risk-neutral valuation framework and as-sume that the expected return is the risk-free rate. The value of the option under the risk-neutral measure is subsequently determined as the discounted expected payoff. The stochastic process underlying the binomial tree can be defined in accordance with the characteristics of the underlying asset, including mean-reversion, drift, stochastic volatil-ity. Another alternative is the use of a Monte Carlo simulation approach. Note that neither binomial trees nor Monte Carlo simulation possesses the analytical traceability.

Since financial markets evolve continuously, the demand for innovation in financial mod-els remains. Tailor-made financial products, growing portfolios and increasing computa-tional power allow quantitative analysts to extend the landscape of mathematical models. Models with the analytical traceability property are of particular interest. Within banks and insurance companies, interest rate term structure models, stochastic volatility op-tion pricing models and credit default models are examples of more recent applicaop-tions of mathematics in finance.

One class of models that is of special interest and has a wide applicability is the so-called affine models. Affine models are a particular class of stochastic models that possesses useful semi-analytical properties and are highly traceable. The main property of interest for financial applications, is that if a stochastic model fulfills the admissability conditions and hence is said to be affine, the (conditional) characteristic function is known. The

(6)

(conditional) characteristic function of a random variable completely describes the (con-ditional) probability distribution. Hence, for pricing financial products, e.g. options, one can use this property. Another advantage is that when affine models are used in simulation studies, e.g. Asset Liability Management (ALM) studies, one does not need to perform nested simulations due to the semi-analytical properties.

The framework to price options using affine models has a generic form. The first step is that one has to determine the underlying financial phenomena to be modelled. There-after, the appropriate underlying process and mathematical form has to be selected. Furthermore, appropriate and sufficient high-quality data is required. Given that the form of the process is affine, the following general set-up is opted for:

• Verify whether the admissibility conditions for the affine form are fulfilled.

• Solve a set of Ordinary Differential Equations (ODEs).

• Define the payoff function for the option.

• Use the solutions of the ODEs to determine the (conditional) characteristic function of the underlying process related to the payoff function.

• Use the Fourier Inversion Theorem to find the probability distribution function.

• Determine the price of the option given the model parameters.

In this thesis, we will explore the class of affine models in light of financial modelling and particularly option pricing. Since the approach possesses complex elements, a more basic model to describe the financial phenomena is considered first, i.e. the Vasicek model. The advantage is that the distribution function is known and hence a closed-form solution for the characteristic function exists. This way, our general framework can be properly explained and validated. Having established a general framework, the models can be extended and more complex structures can be evaluated relatively easily. Although the affine models become more complex, they still have a semi-analytical solu-tion for the characteristic funcsolu-tion. As an example, the framework is extended with the Heston model. The Heston model has an additional stochastic process for the volatility component.

An additional aspect of affine models is parameter estimation. The previous set-up is conditional upon the model parameters. Hence, an extension is to estimate the model parameters using (market) data, e.g. caps and floors, call and put options. By esti-mating the model on current market data, the parameters are determined under the risk-neutral measure. Thereafter, the model could be used for simulation.

The outline of the thesis is as follows. The next chapter discusses the literature used as inspiration for the thesis. In chapter three, the stochastic differential equations are considered. In chapter four, the theory of affine models is introduced. Chapter five links

(7)

affine models to option pricing. In chapter six, the estimation procedure is described. The results are found in chapter seven with the concluding notes in chapter eight.

(8)

2. Literature review

During the research phase different academic papers have been consulted. The initial thoughts related to the final research were found in the papers of Heston (1993), Duffie et al. (2000) and Filipovic (2009).

Research related to option pricing starts with Black and Scholes (1973). They derived a partial differential equation that gives the price of an option for a particular moment in time. The Black-Scholes differential equation can be considered as the most important equation in finance. The no-arbitrage insights of Black and Scholes (1973) formed the building-block for many financial models, e.g. Vasicek (1977), Hull and White (1993), Heston (1993). The no-arbitrage condition relates to the definition of the risk-neutral valuation framework. In this framework one can calculate the value of an asset as the discounted expected future payoff under a particular probability measure, the so-called risk-neutral measure. Here, we also use the risk-neutral valuation framework.

The paper of Heston (1993) has demonstrated us how call option prices can be calculated using numerical inversion techniques applied on the characteristic function. The pricing formula presented by Heston (1993) involves a two-factor model including a process for stochastic volatility. Heston (1993) shows the existence of a closed-form solution for the characteristic function of the model. Moreover, it is demonstrated how the call option price is calculated by inversion of this characteristic function.

The paper of Duffie et al. (2000) extends the foundations of Heston (1993) towards a general framework for affine jump-diffusion processes. Different financial applications are linked to these processes. In the paper, the analytical treatment of the class of Fourier (or Laplace) transforms has been demonstrated. Duffie et al. (2000) shows the existence of a closed-form solution for the characteristic function, in case a model be-longs to the class of affine jump-diffusion processes. As a result, the closed-form solution for the transform results in a semi-analytical pricing expression which can be used for different financial applications.

In Filipovi´c (2009), the affine diffusion processes are revisited and admissability con-ditions are derived. A wide range of applications has been studied. The author presents an extensive study with theoretical insights and financial applications. Amongst others, bond pricing, bond option pricing and stock option pricing are explained making use of different affine models, e.g. Vasicek and Cox Ingersoll Ross (CIR). The theory presented in Filipovi´c (2009) serves as a building block for our research. Moreover, the presented results will be used to verify the framework as described within this thesis.

(9)

3. Stochastic Differential Equations

A stochastic process is defined as an infinite collection of random variables defined on a common probability space. Probability models are defined by three components. First, the sample space, which is a set whose elements are called sample points or outcomes. Second, the class of events with all subsets of the sample space. Third, a probability measure that assigns a nonnegative number to each outcome. In case of throwing a dice, {1, 2, 3, 4, 5, 6} is the sample space, 2, 4, 6 is an even event in the sample space and the corresponding probability measure is 1₂. Stochastic processes are part of Stochastic Differential Equations (SDE). This chapter introduces the SDE and the forms as defined by Vasicek (1977) and Heston (1993).

3.1. General form of Stochastic Differential Equations

The phenomena of the SDE is used to model random processes, e.g. stock prices, volatil-ities, instantaneous interest rates (or short rates). In general a SDE contains a so-called drift term and diffusion term, both depend on the current value of the process only. Moreover, a random white noise term is included in a SDE; which is defined as the derivative of the Brownian motion or Wiener process. The general mathematical repre-sentation of a SDE of X(t) is:

dX(t) = µ(X(t))dt + σ(X(t))dW (t) (3.1) The SDE consists of a drift term µ(X(t)), a Wiener process W (t) and a volatility term σ(X(t)) that describes the diffusion matrix defined by:

Ω(X(t)) = σ(X(t))σ(X(t))> (3.2) The solution of the general SDE form is obtained by integration of both sides of the equation over the interval t to T leading to:

X(T ) = X(t) + Z T t µ(X(s))ds + Z T t σ(X(s))dW (s) (3.3)

For the general solution, one integrates over t to T . One could also integrate over the interval [0, t] or set t = 0 and T = t. In case of a deterministic drift and volatility term, one could in a straightforward way solve the SDE explicitly. In general, the drift and volatility term are stochastic processes which makes finding solutions more challenging.

(10)

process X(T ) enables the calculation of expectations required for pricing purposes. Ex-plicit expressions are rare although expressions in terms of numerical solutions do occur more often.

3.2. The Vasicek Stochastic Differential Equation

As mentioned in the previous section, in some instances SDE models have an explicit analytical solution in terms of the distribution of X. Here, we present an example: the Vasicek model. Other well-known models as Heston, Hull and White or other multi-factor models do not have explicit solutions and hence numerical techniques are required to calculate probabilities and densities.

3.2.1. The distribution of the Vasicek model

The Vasicek model dates back to 1977. The model provides a specification of the term structure of interest rates in the general SDE form with r(t) = X(t), drift term b − βr(t) and constant volatility σ. The mathematical form of this one-factor model is:

dr(t) = (b − βr(t))dt + σdW (t) (3.4) The probability distribution of the interest rate modelled by Vasicek (1977) is known1_.

The Vasicek model contains a stochastic Itˆo integral which is normally distributed with mean zero and a variance that can be determined using the Itˆo isometry theorem. For the Vasicek model, the interest rate under the risk-neutral measure is normally distributed with mean and variance given in equation (3.5). For the proof refer to Appendix A.2.

EQ[r(T )|Ft] = e−βτr(t) + b β 1 − e −βτ V arQ[r(T )|Ft] = σ2 2β(1 − e −2βτ )₎ (3.5)

Using the approach one can also retrieve the distribution of r(0) given the information available at time 0 also referred to as the filtration F0. Note that for notational reason

we have defined τ = T − t.

3.2.2. The (conditional) characteristic function of Vasicek

As outlined in the introduction, the advantage of using affine processes is that the form of the (conditional) characteristic function is known and can be solved numerically. The (conditional) characteristic function of affine processes is described by a set of ODEs. In the next chapter, the formal definitions and admissability conditions related

1

Note that we mention the interest rate, however in many applications the Vasicek model, alike other term structure models, is used to model the short rate. The short rate is the instantaneous spot rate. In other words, the interest rate at which financial institutions lend money for an infinite small period of time.

(11)

to affine processes will be introduced. The conceptual idea is best demonstrated by the (conditional) characteristic function of the Vasicek model. The general form of the (conditional) characteristic function with, fX(·) also conditional, is:

ηX(t) ≡ EeitX|Ft=

Z ∞

−∞

eitXfX(x)dx (3.6)

Since the (conditional) probability distribution for Vasicek is normal with mean µ and variance σ2 _{defined in equations (3.5), the characteristic function is as follows:}

ηX(t) = eiµt− 1 2σ 2_t2 = e e−βτr(t)+βb(1−e −βτ₎_it−σ2(1−e−2βτ )t2 4β (3.7)

In case of the Vasicek model, both the (conditional) characteristic function and the (conditional) distribution function have a known analytical expression. In case one encounters a stochastic model for which only one of the two is known, would it then be possible to derive the other? It results from the Fourier inversion theorem.

3.2.3. The characteristic and probability distribution function

Two distinct distribution functions never have the same characteristic function. The characteristic function uniquely determines the probability distribution. The relation between the characteristic function and the distribution function can be utilized using the Fourier (or Laplace) transform theorem 2. The original application of the Fourier transform was in physics, to transform a function of time into frequency. The definition of the Fourier transform ˆf (ω), dates back to Fourier (1878) and is given by:

ˆ f (λ) =

Z ∞

−∞

e−iλXfX(x)dx (3.8)

The Fourier transform is valid on the complex plane C. The fundamental inversion from the Fourier analysis shows that an inversion formula exists of the form:

fX(x) = 1 2π Z ∞ −∞ eiλXf (λ)dλˆ (3.9)

By choosing t = −λ one can observe that the characteristic function is equal to the Fourier transform: ˆ f (−t) = Z ∞ −∞ eitXfX(x)dx = ηX(t) (3.10)

Therefore, the characteristic function is the Fourier transform with t = −λ. Substituting the characteristic function into the inverse Fourier transform gives us the desired result to retrieve the distribution function. This only holds whenever the characteristic function is known. One can present the inverse Fourier transform as:

2

The Laplace transform is defined as E[esX]. Choosing for s, the term −it in the complex space gives us the Fourier transform.

(12)

fX(x) = 1 2π Z ∞ −∞ e−itXf (−t)dt =ˆ 1 2π Z ∞ −∞ e−itXηX(t)dt (3.11)

One can retrieve the distribution function of a process X using the characteristic func-tion. Note that using the conditional expectations or conditional distribution functions involves conditioning on all information available at time t, the filtration Ft. For the

model of Vasicek and the other models of the affine class, the characteristic function and hence the inverse method will determine the distribution function. For most models a numerical approach for inversion is required but for Vasicek we can demonstrate the derivation of the normal probability density function using the characteristic function analytically. Start with substituting the characteristic function (3.7) in terms of µ and σ2 into the inverse Fourier transform:

fX(x) = 1 2π Z ∞ −∞ e−itXeiµt−12σ 2_t2 dt = 1 2π Z ∞ −∞ e−itX+iµt−12σ 2_t2 dt (3.12)

In order to solve the integral we use the Gaussian integral, which defines that:

Z ∞

−∞

e−a(x−b)2dx =r π

a (3.13)

We define a(t) ≡ it and use the complex characteristic i2 = −1, then:

fX(x) = 1 2π Z ∞ −∞ e−it(X−µ)+12σ 2_(it)2 dt = 1 2π Z ∞ −∞ e−a(t)(X−µ)+12σ 2_a(t)2 dt = 1 2π Z ∞ −∞ eσ22 a(t) 2₋2(X−µ) σ2 a(t) dt = 1 2π Z ∞ −∞ eσ22 a(t)− X−µ σ2 2 −(X−µ)2 2σ2 dt = 1 2πe −(X−µ)2 2σ2 Z ∞ −∞ eσ22 a(t)− X−µ σ2 2 dt (3.14)

The first term outside the integral is part of the normal probability density function. Next, the a(t) terms are replaced by the actual form:

fX(x) = 1 2πe −(X−µ)2 2σ2 Z ∞ −∞ eσ22 it− X−µ σ2 2 dt = 1 2πe −(X−µ)2 2σ2 Z ∞ −∞ e−σ22 t 2_+2itX−µ σ2 − X−µ σ2 2 dt = 1 2πe −(X−µ)2 2σ2 Z ∞ −∞ e−σ22 t 2₋2t i X−µ σ2 − X−µ σ2 2 dt = 1 2πe −(X−µ)2 2σ2 Z ∞ −∞ e−σ22 t− 1 i X−µ σ2 2 dt = 1 2πe −(X−µ)2 2σ2 r 2π σ2 = 1 σ√2πe −(X−µ)2 2σ2 (3.15)

(13)

In the third step, we use that i2 = −1 or i = −1_i. While in the fourth step, we use

1

i2 = 1_i · 1_i = −i · −i = i2 = −1. The solution of the Gaussian integral of (3.13) can

be used by defining x ≡ t, a ≡ σ₂2 and b ≡ 1_iX−µ_σ2 . This way, the probability density

function is derived using the inverse Fourier transform and the characteristic function only. If one substitutes, the mean and variance of the normally distributed Vasicek, see (3.5), the exact distribution can be used for estimation, pricing and simulation purposes.

3.3. The Heston Stochastic Differential Equation

The Vasicek model is a simplified one-factor stochastic model with the underlying as-sumption of constant volatility. To accurately price call and put options in financial markets, stochastic volatility is required. Introducing a separate stochastic process is in accordance with the random nature of volatility, i.e. the dynamics of the volatility is not independent of the level of the underlying. Moreover, stochastic volatility enables to model the skew dynamics of the volatility surface. The Heston model is an example of a two-factor model, where one factor describes the stochastic behavior of the volatility.

3.3.1. The distribution of the Heston model

The Heston model dates back to 1993. The model provides a specification of the stock dynamics in the general SDE form with S(t) = X(t), drift term µS(t) and diffusion term pv(t)S(t). The volatility v(t) follows a separate process. The mathematical form is:

dS(t) = µS(t)dt +pv(t)S(t)dWS(t)

dv(t) = κ(θ − v(t))dt + σv

p

v(t)dWv(t)

(3.16)

The processes dWS(t) and dWv(t) are Wiener processes with correlation ρ. The

stochas-tic process of v(t) is a CIR process. The CIR process has a Noncentral Chi-Squared distribution. The distribution of the Heston SDE, which depends both on S(t) and v(t), is however not explicitly known. Therefore contrary to the Vasicek example, the actual probability distribution function has no familiar form and moreover the characteristic function seems to have no explicit (analytical) solution.

The form of the Heston model may indicate that the distribution of the SDE does not have an explicit form. In the next chapter, the concept of affine models is introduced. Whenever a model satisfies particular admissibility conditions, it is said to be affine. The advantage of affine models is that the characteristic function has an explicit form. With respect to the form of the Heston model (3.16), it turns out that the admissibility conditions are not satisfied. The reason is the S(t) term in the diffusion part of the process. To overcome his, we use a transformed form, i.e. log S(t). The transformed process is found using Itˆo’s lemma. Itˆo’s lemma gives the form of the stochastic process of a function, e.g. log S(t), when the stochastic process of the function variable, e.g. S(t), is known.

(14)

The stochastic process of log S(t) is: d log S(t) = (µ − 1 2v(t))dt + p v(t)ρdW1(t) + p 1 − ρ2_dW 2(t) dv(t) = κ(θ − v(t))dt + σv p v(t)dW1(t) (3.17)

Refer to Appendix A.6 for the proof. Note that we also write the Wiener processes (dWS(t), dWv(t)) with correlation ρ as two independent Wiener processes (dW1(t), dW2(t)).

Due to the transformation, the diffusion part of the log S(t) process only depends on the stochastic volatility v(t) component and not on S(t) anymore.

(15)

4. Properties of affine models

In the previous chapter the relation between the characteristic function and the dis-tribution function was explained by example of Vasicek. Using the Fourier inversion theorem, one can relate the probability distribution function to the characteristic func-tion. For affine models, the characteristic function has an explicit form. In this chapter, the admissibility conditions and properties of affine models are considered. Moreover, a numerical method to solve the Riccati equations is proposed when analytical solutions lack. The method is demonstrated by example of the Vasicek and Heston model.

4.1. Formal definition of the affine processes

In Filipovi´c (2009), the class of affine processes is described. The application of pricing financial products like options using affine models is introduced in Duffie and Singleton (2000) where the initial set-up was provided by Heston (1993). Following the theorems of Filipovi´c (2009), a SDE X(T ) is called affine if the characteristic function conditional upon the filtration at time t, is of an exponential affine form1 in X(t):

eu>X(T )

Ft= eφ(T −t,u)+ψ(T −t,u)

>_X(t)

(4.1) Here φ(T − t, u) and ψ(T − t, u) are functions in the complex space and solutions of a system of ODEs, named Riccati equations. The equations hold for all u ∈ iRd; where d is the dimension of the SDE. Whenever X(T ) is affine, the diffusion matrix and drift are also affine. The argumentation is invertible. For the proof refer to Filipovi´c (2009). Suppose the diffusion matrix and drift of SDE X(T ) are affine and a solution exists for the system of ODEs such that φ(t, u) + ψ(t, u)TX(t) has a non-positive real part for all t ≤ 0 and u ∈ iRd, then X(T ) is affine with (conditional) characteristic function of the form (4.1)

A SDE on the canonical state space is affine if and only if the parameters of the diffusion matrix and drift vector fulfill the admissability conditions and there exists a solution of the system of Riccati equations2. The formal definition is given in Filipovi´c (2009):

Theorem 1: The process X on the canonical state space Rm

+ × Rn is affine if and

only if the drift and the diffusion of the process are affine, i.e. Ω(X(t)) = a +Pd

i=1Xiαi

and µ(X(t)) = b +Pd

i=1Xiβi for a, αi, b, and βi which are admissible in the following

1_{An affine form is characterized as a + bX} 2

Riccati equations are first-order ordinary differential equations which have a quadratic nature. The general form is ∂tH(t, u) = AH(t, u)2_{+ BH(t, u) − C.}

(16)

sense:

1. a, αi symmetric and positive semi-definite

2. aII = 0 and thus aIJ = a>J I = 0

3. αj = 0 for all j ∈ J

4. αi,kl= αi,lk = 0 for k ∈ I \{i}, for 1 ≤ i, l ≤ d

5. b ∈ Rm+ × Rn

6. BIJ = 0

7. BII has positive off-diagonal elements

In this case, the corresponding system of Riccati equations simplifies to: ∂φ(t, u) ∂t = 1 2ψJ(t, u) >_a J JψJ(t, u) + b>ψ(t, u) φ(0, u) = 0 ∂ψi(t, u) ∂t = 1 2ψ(t, u) > αiψ(t, u) + βi>ψ(t, u) ∂ψJ(t, u) ∂t = B > J JψJ(t, u) ψ(0, u) = u (4.2)

and has a unique global solution (φ(·, u), ψ(·, u)) for all initial values of u. The index sets I and J are defined as I = {1, ..., m} and J = {1, ..., n}. Moreover a and αi are

d × d-matrices and b and βi are vectors of length d × 1. The dimension d is the sum of

the dimension of the state variables n and the stochastic volatility components m, i.e. d = m + n. The parameters a, αi, b, and βi are such that X does not leave the set (or

state space) X where X = Rm₊ × Rn_.

One can thus verify if a process is affine by writing the process in the general form in accordance with equation (3.1) and subsequently use Theorem 1 to assess the parameters a, αi, b, and βi.

4.2. A solution for the (conditional) characteristic function

The (conditional) characteristic function of an affine process X is completely described by the solutions of the system of Riccati equations and the information at time t. As we

(17)

will demonstrate in the next section, the Riccati equations of the Vasicek model have an analytical form. However, for most processes with an extended form, e.g. multiple factors, stochastic volatility, the system of equations can only be solved numerically. The numerical approach is explained hereafter. The validity of the approach is demonstrated and validated by the Vasicek example.

4.2.1. An iterative method for solving the Riccati equations

Given the affine representation of the process as in equation (3.1) and suppose the admis-sibility conditions are satisfied, the Riccati equations are of the form (4.2). The system of differential equations of Riccati can be solved using the power series of the functions ψi(t, u) and φ(t, u). The power series allows us to find an iterative expression for the

Riccati equations. The proof is given in Appendix A.5.

Theorem 2: The system of Riccati equations of the form (4.2) can be solved nu-merically using the power series:

ψi(t, u) = ∞

X

k=0

Ck,i(u)tk= C0,i(u) + C1,i(u)t + C2,i(u)t2+ ...

φ(t, u) =

∞

X

k=0

Dk(u)tk= D0(u) + D1(u)t + D2(u)t2+ ...

(4.3)

The corresponding numerical solution for ψi(t, u) and φ(t, u) then reads:

Ck+1,i(u) = 1 2u > αiu + βi>u (k = 1) Ck+1,i(u) = 1 k + 1 1 2 k X j=0

C_j>(u)αiCk−j(u) + βi>Ck(u)

(k > 1) Dk+1(u) = 1 2u >_{au + b}>_u _{(k = 1)} Dk+1(u) = 1 k + 1 1 2 k X j=0

C_j>(u)aCk−j(u) + b>Ck(u)

(k > 1)

(4.4)

Note that a and αi are square matrices and b and βi are vectors with the dimensions

depending on the number of factors in the affine process.

4.2.2. Iteratively solving the (conditional) characteristic function

The numerical method to solve the Riccati equations expresses the solution for ψi(t, u)

and φ(t, u) as power series. An approximation error is introduced during the implemen-tation caused by the infinite power series being approximated by a finite power series with boundary N < ∞. Remark that the power series depends on the time t. Suppose that one would like to price a financial instrument using the (conditional) characteristic

(18)

function of an affine process, see equation (4.1). In case the maturity of the instrument and hence the time component in the power series is relatively large, the approximation error of the numerical method will also be relatively large. Moreover, increasing the maturity, ceteris paribus, increases the approximation error.

The exponentially affine form of the (conditional) characteristic function allows us to as well solve this expression in an iterative way. If one returns to the original result of the process being affine in equation (4.1), one observes that for u ∈ iRdand standing at time t = 0, the (conditional) characteristic function for our process X is:

E h eu>X(T ) F0 i = eφ(T,u)+ψ(T,u)>X(0) (4.5) The maturity T can be divided into M sub-intervals of equal size ∆. Then the Law of Iterated Expectations, which tells us that E[X|I1] = E[E[X|I2]|I1], allows us to write

equation (4.5) as: Eheu>X(T ) F0 i = EhEheu>X(T ) FT −∆ i F0 i (4.6)

Using the affine (conditional) characteristic function with t = T − ∆ and the fact that T − (T − ∆) = ∆ gives: Eheu>X(T ) F0 i = Eheφ(∆,u)+ψ(∆,u)>X(T −∆)|F₀i (4.7)

Bringing the constant term eφ(∆,u) outside the expectation and rename u(1) ≡ ψ(∆, u) results in a familiar form:

Eheu>X(T ) F0 i = eφ(∆,u)Eheu(1)>X(T −∆)|F₀i (4.8)

Repeating the step with T = T − ∆ and u(2)≡ ψ(∆, u(1)_{) results in:}

Eheu>X(T ) F0

i

= eφ(∆,u)Eheφ(∆,u(1))+ψ(∆,u(1))>X(T −2∆)|F₀i = eφ(∆,u)+φ(∆,u(1))E

h

eu(2)>X(T −2∆)|F0

i (4.9)

And for all M steps on the interval the resulting expression will be:

E h eu>Y (T ) F0 i = eφ(∆,u)+φ(∆,u(1))+...+φ(∆,u(M ))E h eu(M )>Y (0)|F0 i (4.10)

The number of intervals M , hence the size of ∆, has to be chosen in an appropriate way to use the computation time efficiently and at the same time increase the chance of finding a solution. Note that after each iteration the u variable which is the start value of the iteration is updated, by the next start value, i.e. ψ(∆, u), ψ(∆, u(1)), ..., ψ(∆, u(M −1)). The last iteration states that u(M ) ≡ ψ(∆, u(M −1)_).

(19)

4.3. The numerical method in practice with examples

In the previous subsections, the theoretical framework for numerically solving the Riccati equations has been presented. The numerical method has been implemented and tested in MatLab. The solver follows the iterative approach. In other words, suppose the Riccati equations have to be solved for a contract with maturity T , one divides the interval in M subintervals of equal size ∆ and solves the Riccati equations accordingly. Note that for the first interval [T − ∆, T ], the starting value of the method is u (see equation (A.27). For the next interval, the starting value is updated to ψ(∆, u) and T = T − ∆. Iterating over all intervals results in the ultimate numerical solution.

4.3.1. The solution of the Riccati equations for Vasicek

In order to validate the Riccati solver, the closed-form solutions for the Riccati equations involving the Vasicek model have been implemented. The form of the Vasicek model, equation (3.4), can be written in the general (affine) SDE form (3.1):

µ(r(t)) = b − βr(t)

Ω(r(t)) = σ2 (4.11)

The affine representation (4.1) reads as X(t) = r(t), a = σ2, α = 0, b = b and β = −β. The parameters fulfill the admissibility conditions.

Figure 4.1.: The graphs show the difference between the numerical and analytical method for solving the Riccati equations of the Vasicek model. The left graph shows the abso-lute differences related to φ(t, u) and the right graph shows the absoabso-lute differences related to ψ(t, u) for 100 random u numbers in the complex space. The dimension of the power series is N = 13. Remark the the size of the difference.

Using Theorem 1, the (conditional) characteristic function is given as an exponential expression of the solution of the Riccati equations. The system of Riccati equations for the Vasicek model is given in the Appendix A.3. Moreover, the analytical solutions are

(20)

derived in the Appendix. Note that the dimension of the SDE is d = 1 and a single solution for ψi(t, u) = ψ(t, u) occurs. Vasicek assumes a constant volatility and hence

no separate stochastic volatility process is defined. The solutions are:

ψ(t, u)) = ue−βt φ(t, u) = 1 4β σ2u2 1 − e−2βt + 4bu 1 − e−βt (4.12)

The Riccati solutions φ(t, u) and ψ(t, u) are considered using our numerical solver and the closed-form expressions. The absolute difference is calculated for both methods using a random set of complex u values. For the numerical method, the infinite summation of the power series has to be bounded. The boundary used is N = 13. In other words, the orders of the polynomial beyond 13 are ignored. The results show an error in the order of 10−13. The relatively small differences confirm the validity of the numerical approach and implementation. Note that increasing the bound, N > 13, further reduces the error however at the same time increases the computation time.

4.3.2. The solution of the Riccati equations for Heston

For affine models more complex than the Vacicek model, the underlying distribution can often not be derived analytically. Moreover, the system of Riccati equations can only be solved using numerical methods. For the transformed Heston model (3.17), we verify the admissability conditions. Note that from Theorem 1, one observes that the dimensions of I and J are both one, hence d = 2. The index set I relates to the v(t) process 3 and J relates to the log S(t). Denote the variables as:

X(t) =X1(t) X2(t) =log S(t) v(t) (4.13)

Write the model in the general form (3.1) with µ(X(t)) and Ω(X(t)) in terms of a, α1,

α2 b, β1 and β2 originating from Theorem 1.

b = µ κθ β1 = 0 0 β2 = −1 2 −κ a = 0 0 0 0 α1 = 0 0 0 0 α2 = 1 ρσv ρσv σv2

The (affine) admissability conditions are fulfilled, see Appendix A.3 for the references. Note that the matrix α1 and the vector β1 are filled with zeros, hence the system of

Riccati equations reduces to finding the solution of φ(t, u) and ψ2(t, u). The solution

of ψ1(t, u) is according to the system of Riccati equations (4.2) given by ∂ψ1_∂t(t,u) = 0.

Together with the boundary condition, the solution reads ψ1(t, u) = ψ1(0, u) = u1. The

Riccati equations for φ(t, u) and ψ2(t, u) are respectively:

(21)

∂φ(t, u) ∂t = 1 2ψ2(t, u) · 0 · ψ2(t, u) +µ κθ ψ1(t, u) ψ2(t, u) ∂ψ2(t, u) ∂t = 1 2ψ1(t, u) ψ2(t, u) 1 ρσv ρσv σv2 ψ1(t, u) ψ2(t, u) +−1 2 −κ ψ1(t, u) ψ2(t, u) (4.14)

The system of Riccati equations can be solved using the numerical approach. However, one could also derive an analytical solution. The derivation is presented in Appendix (A.7) and uses the Lemma 10.12 of Filipovi´c (2009).

ψ2(t, u) = −

u1(1 − u1)(eθt− 1) − (θ(eθt+ 1) + (ρσvu1− κ)(eθt− 1))u2

θ(eθt_{+ 1) − (ρσ} vu1− κ)(eθt− 1) − σv2(eθt− 1)u2 (4.15) φ(t, u) = µu1t + 2κθ σ2 v log 2θe θ−(ρσv u1−κ) 2 t θ(eθt_{+ 1) − (ρσ} vu1− κ)(eθt− 1) − σv2(eθt− 1)u2 (4.16)

In the analytical solution of the Riccati equations of Heston, the θ is a function of the parameters (σv, ρ, κ, u1) and has the familiar form:

θ =p(ρσvu1− κ)2+ σv2u1(1 − u1) (4.17)

Also for the Heston model, the Riccati solutions for φ(t, u), ψ1(t, u) and ψ2(t, u) are

calculated using the numerical solver and the analytical solution. The absolute difference is calculated for both methods using a random set of complex random u values. The results for φ(t, u), ψ1(t, u) and ψ2(t, u) show relatively small errors in the order of 10−9.

Figure 4.2.: The graphs show the difference between the numerical and analytical method for solving the Riccati equations of the Heston model. The left graph shows the absolute differences related to φ(t, u), the middle graph shows the absolute differences related to ψ1(t, u) and the right graph shows the absolute differences related to ψ2(t, u) for

100 random u numbers in the complex space. The dimension of the power series is N = 13. Remark the size of the difference.

One interesting feature has been encountered when comparing the numerical and ana-lytical solutions related to φ(t, u). In case the complex numbers u increase in size, e.g.

(22)

1000 + 1000i, the solution is not unique. The imaginary part of the solution is a multiple of 2πi. For the calculations and option pricing, the non-uniqueness has no impact, since e2πi = 1, e.g. the characteristic function provides the same results whenever a complex number is used or a multiple of 2πi.

4.4. Discounting incorporated in the affine models

One particular financial variable also modelled is the short rate or instantaneous spot rate. The class of affine models is well-suited for modelling the short rate. Given the general form of the affine process X(t), one can define the short-rate process as:

r(t) = c + γ>X(t) (4.18) Where the parameters c and γ are respectively defined in R and Rd. Note that by choosing c = 0, γ = 1, the Vasicek model can be used for modelling the short rate. The extension towards the short-rate process results in a highly useful property in accordance with the (conditional) characteristic function earlier. In Filipovi´c (2009) it has been proven that the conditional risk-neutral expectation of the discounted version of the function euTX(T ) has an affine representation:

EeRtTr(s)dseu >_{X(T )}

Ft= eΦ(T −t,u)+Ψ(T −t,u)

>_X(t)

(4.19) Where Φ(t, u) and Ψ(t, u) are solutions of a system of Riccati equations:

∂Φ(t, u) ∂t = 1 2ΨJ(t, u) > aJ JΨJ(t, u) + b>Ψ(t, u) − c ∂Ψi(t, u) ∂t = 1 2Ψ(t, u) >_α iΨ(t, u) + βi>Ψ(t, u) − γi ∂ΨJ(t, u) ∂t = B > J JΨJ(t, u) − γJ (4.20)

Also here ψ(0, u) = u and φ(0, u) = 0 have to be satisfied. Note the additional terms c and γ in the system of Riccati equations. Moreover, a different notation is used Φ(t, u) and Ψ(t, u) instead of φ(t, u) and ψ(t, u) for the processes without discounting in (4.2). In Appendix A.4, we derived the analytical solutions of the Riccati equations for the Vasicek model including discounting. The solutions are more extensive compared to equations (4.12): Ψ(t, u) = ue−βt+e −βt_{− 1} β (4.21) Φ(t, u) = 1 2σ 2u2 2β(1 − e −2βt_{) +} u β2(1 − e −2βt_{) −}2u β2(1 − e −βt_{) +} 1 2β3(1 − e −2βt₎ − 2 β3(1 − e −βt_{) +} t β2 +bu β(1 − e −βt_{) +} 1 β2(1 − e −βt_{) −} t β (4.22)

(23)

Figure 4.3.: The graphs show the difference between the numerical and analytical method for solving the Riccati equations of the Vasicek model for short rate modelling hence including discounting. The left graph shows the absolute differences related to φ(t, u) and the right graph shows the absolute differences related to ψ(t, u) for 100 random u numbers in the complex space.

Comparing the analytical solutions of the Riccati equations with our numerical method for 100 random complex numbers show differences of the order no larger than 10−8.

With the solutions of the Riccati equations and the affine representation (4.19), one can easily determine the price of a discount bond with maturity T . The price of the bond results whenever u = 0 such that eu>X(T )_{= e}0 _{= 1 and hence:}

P (t, T ) = EeRtTr(s)dsF_t= eΦ(T −t,0)+Ψ(T −t,0) >_X(t)

(4.23) Note that P (t, T ) represents the price of a bond with maturity T at time t. Occasionally, one will encounter that Φ(T − t, 0) and Ψ(T − t, 0) are defined as new variables. For example, A(T − t) ≡ Φ(T − t, 0) and B(T − t) ≡ Ψ(T − t, 0). Under this definition, the bond price reads as P (t, T ) = e−A(T −t)−B(T −t)>X(t).

(24)

5. Pricing options using affine models

Since the (conditional) characteristic function of an affine SDE has an explicit expres-sion in terms of the Riccati equations, see equation (4.1), one can use the technique of Fourier inversion1 to retrieve the probability distribution and from there the cumulative distribution. Having the cumulative distribution function of the process underlying the financial product, one can determine the exact value. In case of an option, one can determine the probability that the option expires in-the-money. In this chapter, the mathematical building blocks for option pricing are considered. The building blocks combined enables us to price both European call and put options.

5.1. The Fourier inversion method

The characteristic function for affine models is known by definition. Since the character-istic function uniquely defines the probability distribution function, one is able to retrieve this function and use it for calculating expectations, for example in case of option pricing. Different techniques are available for numerically inverting the characteristic function. In this thesis, the numerical inversion method introduced in Den Iseger (2006, 2008) is used. This method uses a well-known technique for the inversion of the transform. The transform is in fact the characteristic function.

The inversion method relates the values of the probability distribution function to the Laplace (or Fourier) transform series using the Poisson Summation Formula. The nu-merical inversion method has been extended in the later paper Den Iseger (2008), such that the method returns the function values for an arbitrary grid of points instead of a uniform grid (or single point). Den Iseger (2008) first establishes a relation between the Laplace (Fourier) series in the coefficients of the Legendre expansion2 _{and the}

orig-inal Laplace (Fourier) transform series. The infinite Laplace (Fourier) transform series is approximated by a finite series using the Gaussian quadrature rule. Having a finite series, one can determine the coefficients of the Legendre expansion using the Inverse Fast Fourier Transform (IFFT). The coefficients are subsequently used to determine the function values of the probability distribution function.

Note that for affine models the values of the Laplace (Fourier) transform series are

1

The Fourier inversion theorem states that a well-behaved function can be recovered from its Fourier transform.

2_{The Legendre polynomials are a mathematical technique to make a polynomial approximation for a}

certain function on a particular interval. Another often encountered polynomial approximation is the Taylor approximation.

(25)

available since the characteristic function is known. Recall the similarity between the Fourier transform and the characteristic function. This allows us to determine the co-efficients of the Legendre expansion using IFFT. Having calculated the coco-efficients, one can determine the function values of the probability distribution function. For the math-ematical details, which are beyond the scope of the thesis, we refer to the Den Iseger (2006) and (2008) papers.

5.2. From the CCF to CDF

The Fourier inversion of the (conditional) characteristic function of an affine process results in the probability distribution of process X. Since one needs the cumulative dis-tribution in order to calculate the expectation under the risk-neutral probability mea-sure, the technique of convolution is used. The definition of the cumulative distribution function (CDF), P [X ≤ x], is explicitly expressed as:

F (x) = E[1{X≤x}] = Z ∞ −∞ 1{X≤x}f (x)dx = Z ∞ −∞ 1{x−X≥0}f (x)dx (5.1)

The theory of convolution states that the convolution of two functions, denoted by ∗, is the integral of the product of the two functions after one is reversed and shifted. One can define the indicator function in equation (5.1) as the shifted function g(x − X) and the probability density function as f (X). The convolution then implies:

Z ∞

−∞

g(x − X)f (X)dx = (f ∗ g)(x) = f (x) ∗ g(x) (5.2)

The Fourier inversion gives us the probability density function f (x). The probability distribution g(x) can explicitly be calculated by, for example, the Laplace inversion:

ˆ L(ω) = Z ∞ 0 e−ωxg(x)dx = Z ∞ 0 e−ωx1{x≥0}dx = Z ∞ 0 e−ωxdx = 1 ω (5.3) Concluding in case of the affine models, the inversion technique results in the (condi-tional) probability distribution function f (x) (for notational convenience, we omit the conditioning on Ft). In order to find the cumulative distribution, one needs to multiply

the (conditional) probability distribution function by the result of equation (5.3), which is the same as dividing by ω. When using the Fourier inversion theorem, ω are the (complex) inversion numbers.

5.3. Change of numeraire

Under the risk-neutral valuation framework, one calculates the value of an asset as the discounted expected future payoff under a particular probability measure, the so-called risk-neutral measure Q. Under the risk-neutral probability measure, the expectations are relatively to the money market account. This relative measurement or benchmark

(26)

is often referred to as the numeraire. The general representation of the arbitrage-free price of a claim ξ(X(T )) under the risk-neutral measure is:

π(t, T, K) = Ee−RT

t r(s)ds_{ξ(X(T ))}F_t (5.4)

Instead of using the risk-neutral measure and hence the money market account as nu-meraire, one can also use another tradeable asset as numeraire. In case of option pricing models for example, the discount bond with maturity T is often used as numeraire. This numeraire is referred to as the T -forward measure. The arbitrage-free price under the T -forward measure is defined as:

π(t, T, K) = P (t, T )E_QTξ(X(T )) (5.5)

The theory to use when one would like to change the numeraire is the so-called Radon-Nikodym derivative. The Radon-Radon-Nikodym derivative allows one to change the numeraire. For example, the expectation under the risk-neutral measure can be written as an expec-tation under T -forward measure and vice versa using the Randon-Nikodym derivative. The derivative is an expression of the pricing functions related to the money market account and the zero-coupon bond with maturity T and time t.

5.4. The arbitrage-free option value function

Now that the distribution of process X can be derived using Fourier inversion techniques jointly with convolution, one can observe this in light of option pricing. Consider a financial claim defined as ξ(X(T )). Note that the definition indicates that the claim depends on the state X at time T . The general expression of the arbitrage-free price at time t ≤ T with respect to claim ξ(X(T )) is given in equation (5.4). Note that the most straightforward claim corresponds to a discount bond where ξ(X(T )) = 1. The claim that is considered in light of option pricing is that of a call option on a stock ξ(X(T )) = (S(T ) − K)+. The arbitrage-free price of a call option on a stock at time t, becomes: π(t, T, K) = Ee−RT t r(s)ds(S(T ) − K)+F_t = Ee−RT t r(s)ds₁_{{S(T )≥K}}_{(S(T ) − K)}F_t = Ee−RT t r(s)ds1_{{S(T )≥K}}S(T )F_t−KEe− RT t r(s)ds1_{{S(T )≥K}}F_t (5.6)

The first part of the pricing function depends on more than one stochastic component. In order to find a useful expression, one can use the so-called Esscher transform which is often used in actuarial science. The Esscher transform uses a new probability measure to define the expectation towards this new probability measure in terms of the risk-neutral probability measure, i.e. _E[eehXhX_]. The random variable X has an Esscher transform with

(27)

fX(x; h) = ehXfX(x) E[ehX_] = ehXfX(x) MX(h) (5.7)

Where MX(h) is the moment generating function defined as

R∞

−∞eh

T_x

fX(x)dx. Observe

the similarity between the moment generating function and the characteristic function from equation (3.6) whenever MX(iu) = ηX(u). Writing the Esscher transform of

equa-tion (5.7) in terms of expectaequa-tions results in:

EA[fX(x; h)] =

E[ehXfX(x)]

E[ehX_] (5.8)

The expectation on the left side is with respect to the Esscher probability measure. With respect to the denominator on the right side, recall the form of the Fourier and Laplace transforms from previous sections and the resemblance with the moment generating function.

5.5. The actuarial Esscher transform illustrated

In order to show the application of the Esscher transform, assume that one has two candidate models to price call options: the Vasicek and the Heston model. The Vasicek model assumes constant volatility, while the Heston model has a stochastic volatility component. Hence, the Heston model describes the price of a stock using stochastic volatility. In case of Vasicek, the price of the stock is defined as S(T ) = eX(T ) where X(t) follows the Vasicek process (3.4) 3. In case of Heston, the price of the stock is defined as S(T ) = eX1(T )_{. Note that the Heston model is described by two factors}

(X1, X2). Although the price of the stock only depends on state X1, state X1 depends

on state of X2. Hereafter, the expression for the option price for Heston is shown as

example. Remark the equivalence in case of the Vasicek model with X1(T ) replaced by

X(T ). The arbitrage-free price of the call option using the Heston model becomes in accordance with equation (5.6):

π(t, T, K) = e−r(T −t) E1_{e_{X1(T )}_≥K}eX1(T )_F t−KE1{eX1(T )≥K} Ft = e−r(T −t)E1{X1(T )≥log(K)}e X1(T )_F t−KE1{X1(T )≥log(K)} F_t (5.9)

Note that for option pricing, constant interest rates r(t) = r are assumed. One could extend the framework by introducing a separate stochastic process for interest rates.

As shown before, the Heston model is affine with X(t) = (X1(t), X2(t)). One can easily

observe that by choosing u = [u1, 0], it follows that eu

>_{X(T )}

= e[u10][X1(T ) X2(T )]> ₌

eu1X1(T ) _{= e}u1log S(T )_{. Since the stock price is defined as e}X1(T )_{, the affine form with} 3

Note that we use X(t) instead of r(t) since it involves a process for the stock price and moreover r(t), the interest rate, is assumed fixed for call option pricing, see equation (5.9).

(28)

u = [u1, 0] can be used to determine the characteristic function and thereby the

distri-bution function and expectations related to the X1 process.

The second term of the arbitrage-free call option price in equation (5.9) can be cal-culated using the CDF. The CDF results from inverting the (conditional) characteristic function for X1 with u = [u1, 0]. The first term is more demanding since the expectation

involves two stochastic components. Here, the Esscher transform is used. Recall the form of the Esscher transform and recognize the similarity between the terms in the call option price EeX1(T )₁

{X1(T )≥log K}

Ft and in the Esscher transform E[ehXfX(x)|Ft].

The explicit form of the transform becomes:

EA[fX1(x; 1)|Ft] =

E[eX1_f

X1(x)|Ft]

E[eX1|F_t] (5.10)

Using the Esscher transform, one is able to split the first term of the arbitrage-free call option price in the product of two components EA[fX1(x; 1)|Ft] and E[e

X1_|F

t]. The nice

property of the expectation under the Esscher probability measure is that it can also be expressed in terms of the risk-neutral probability measure, see for example Gerber and Shiu (1994), EA[fX1(x; 1)|Ft] = E[fX1(x; 2)|Ft]. This allows us to calculate the option

(29)

6. Estimation method

The Esscher transform combined with equation (5.9) and the Fourier inversion technique, to retrieve the probability distribution of the underlying process, allows us to price call options. If one would like to price a put option, the put-call parity can be used. The (conditional) characteristic function of an affine model is known and therefore we have the tools to determine the distribution function and hence expectations. In this chapter, the estimation framework is outlined. We start by validating our framework using results from Filipovi´c (2009). Thereafter, the estimation method is discussed.

6.1. The validity of our implementation

Filipović (2009) serves as building block for our thesis and moreover provides numerical examples. Having implemented the pricing framework for affine models in MatLab, we first use the Heston parameter set-up of Filipović (2009) and calculate the options prices. The option prices in Filipović (2009) are compared with our results. This way we have, besides the comparison of our numerical Riccati solving procedure with the analytical solutions, an additional validation of our framework.1

T \ K 0.80 0.90 1.00 1.10 1.20 0.20 0.2016 0.1049 0.0348 0.0074 0.0012 0.40 0.2037 0.1120 0.0478 0.0168 0.0053 0.60 0.2061 0.1183 0.0571 0.0245 0.0100 0.80 0.2088 0.1239 0.0646 0.0311 0.0144 1.00 0.2115 0.1291 0.0711 0.0368 0.0186 T \ K 0.80 0.90 1.00 1.10 1.20 0.20 0.2016 0.1049 0.0348 0.0074 0.0012 0.40 0.2037 0.1120 0.0478 0.0168 0.0053 0.60 0.2061 0.1183 0.0571 0.0245 0.0100 0.80 0.2088 0.1239 0.0646 0.0310 0.0144 1.00 0.2115 0.1291 0.0711 0.0368 0.0186

Table 6.1.: The table shows European call option prices for different strikes and maturities. The Heston model parameter values of Filipovi´c (2009) are used. The notation of the thesis is used where (µ, κ, θ, ρ, σv) = (0.01, 2.00, 0.01, 0.5, 0.10) and the values of

the state variables are equal to (log S(0), v(0)) = (0.00, 0.04). The upper table are the thesis results, the lower table corresponds with Filipovi´c (2009).

1

Note that Filipovi´c (2009) uses numerical integration to price options, where we use the Esscher transform, equation (5.9) and the Fourier inversion method based on Den Iseger (2006, 2008).

(30)

The option prices we calculate, are European call options. The formula is derived from equation (5.4) with t = 0 and a constant interest rate of r = 0.01. The formula reads as:

π(T, K) = e−rTE[(S(T ) − K)+|F0] (6.1)

The Filipović (2009) table of interest is the one with call option prices, i.e. Table 10.2. Using the same parameter set-up, the European call option prices are calculated2. The results are depicted in the Table 6.1. Alike Filipović (2009), we depict the results in four decimals. One can observe that the results match the results of Filipović (2009) closely. For the call option with maturity 0.80 and strike 1.10 the difference in price is 0.0001, while for the other call options the difference in price is 0.0000. The results provide us with an additional verification of our implementation.

Although for the Vasicek model no European call option prices are given in Filipovi´c (2009), we verify the implementation by choosing a set of parameter values for (a, b, β), i.e. (0.04, −0.20, 0.01). Using these parameters, the European call option prices are calculated and depicted in Table 6.2. The option prices are deemed reasonable.

T \ K 0.80 0.90 1.00 1.10 1.20 0.20 0.2055 0.1097 0.0380 0.0074 0.0008 0.40 0.2116 0.1213 0.0544 0.0185 0.0048 0.60 0.2179 0.1315 0.0669 0.0284 0.0102 0.80 0.2241 0.1405 0.0772 0.0372 0.0157 1.00 0.2300 0.1484 0.0861 0.0450 0.0214

Table 6.2.: The table shows European call option prices for different strikes and maturities. The Vasicek model parameter values are used defined as (a, β, b) = (0.04, −0.20, 0.01) and the values of the state variables equal to log S(0) = 0.00.

In this section, European call option prices have been calculated given the model pa-rameters. For applications in financial markets, the opposite approach is more valuable. Given the market prices, one would like to estimate model parameters. Subsequently, one can use the parameters to price other products or conduct Monte Carlo simulations.

6.2. Estimation of the affine models

6.2.1. The set-up and the market data

The next step is to estimate the model parameters in our framework. We start with the Vasicek model and assume that the stock price is described by S(T ) = eX(T )_{. Hence}

X(t), which is then the logarithm of the stock price, follows the Vasicek process (3.4). The solution of the Riccati equations (4.12) gives us the (conditional) characteristic function. The probability distribution function follows from the Fourier inversion. The

2

The definition of the Heston model in Filipovi´c (2009) differs from ours, see Appendix A.6 for the parameter mapping.

(31)

option price π(t) is determined by (5.6). Note that π(t) is a function of the strike K, the maturity T , the interest rate r and the model parameters. Therefore, we will express the option price as π(T, K, a, α, b, β). In case of Vasicek, we have α = 0. The parameter t is omitted since we only consider prices at t = 0.

We downloaded AEX call option data for the month of December 2015. We will de-note the market price of a call option with strike K and maturity T as C(T, K). The strikes of the call options range from 200, 240, 280, 300, 320, 350, 360, 400, 420, 440, 450, 460, 480, 500, 550, 600. The call option maturities are June 2016, September 2016, December 2016, June 2017, December 2017, December 2018 and December 2019. The affine model parameters are chosen such that the sum of squares between market and model prices is minimized.

We use the Levenberg-Marquardt algorithm3 for non-linear functions to optimize the affine parameters a, α, b and β given a set of AEX option prices. In other words, the algorithm minimizes the squared difference between market and model prices given the affine parameters. The optimization routine sums over all strikes and maturities. Note that the objective function is defined in terms of the affine parameters. The results are easily transformed in terms of the actual model parameters, e.g. for the Vasicek model one uses σ =√a, α = 0, b = b and β = −β.

L = min a,α,b,β N X i=1 N X j=1 C(Ti, Kj) − π(Ti, Kj, a, α, b, β) 2 (6.2)

In the estimation routines it is important to maximize the chance of finding the global optimum. Whenever the number of observations is too low or whenever the number of model parameters is too high, the chance of finding a local optimum is relatively high.

6.2.2. The optimal parameter values

The estimation routine requires starting values for the parameters, i.e. a0, b0, β0. In

case the estimation set and hence the number of observations is relatively low, different starting values may result in different optimal parameter sets. In other words, the routine may converge to different local optima. One can verify if the global optimum is found by conducting the estimation using different starting values.

Vasicek Model Set 1 Set 2

Starting values (0.0100, 0.5000, 0.0200) (0.4000, -1.1000, 5.0000) Optimal parameters (0.0338, -0.0516, 0.2838) (0.5149, -2.9391, 17.6481)

Sum of squares 3.1646 3.1646

Table 6.3.: The December 9, 2015 estimation results for the December 2018 call options for seven strikes. For the estimation two sets of starting values for a0, b0, β0are used.

(32)

To illustrate the occurrence of different optima, suppose we observe the market prices of December 9, 2015 for the call options with expiry on December 2018. We use the market prices for different strikes. The objective function (6.2) is used to minimize the sum of squares given the affine parameters. Two different sets of starting val-ues in the estimation routine are used, (a(1)₀ , b(1)₀ , β₀(1)) = (0.0100, 0.5000, 0.0200) and (a(2)₀ , b(2)₀ , β₀(2)) = (0.4000, −1.1000, 5.0000). The optimal parameters and sum of squares are presented in Table 6.3.

The estimated model and market prices for the AEX call options per strike are depicted in Table 6.4. The results are deemed good, even though Vasicek assumes constant volatil-ity. The reason is that for this estimation only options with the same maturity are used. Note that different starting values have led to different optimal parameter estimates, while the estimated prices are practically the same. This indicates the non-uniqueness of the optimum caused by the low number of observations.

Label Strike Market price Price set 1 Price set 2 AEX 12/18 C300 300 124.9200 125.8701 125.8700 AEX 12/18 C350 350 89.1800 88.4879 88.4879 AEX 12/18 C400 400 60.1900 59.2088 59.2088 AEX 12/18 C450 450 38.1200 37.9915 37.9915 AEX 12/18 C500 500 23.1200 23.5672 23.5672 AEX 12/18 C550 550 13.6100 14.2399 14.2399 AEX 12/18 C600 600 7.9800 8.4349 8.4349

Table 6.4.: The table shows the estimation results for the December 2018 call options at De-cember 9, 2015. The estimation is performed with two sets of starting values for a0,

b0, α0. The first column gives the label, the second column the strike with respect

to the AEX index, the third column are the market prices while the fourth and fifth column show the estimated prices.

As the results of Table 6.3 demonstrate, different sets of optimal parameters are found depending on the starting values of the estimation. Hence it is important to verify the estimation results by, for example, using different starting values. Moreover, sufficient market prices are required to find a solution that converges to a global optimum.

In the next chapter, we will estimate the model parameters based on the full set of observed market prices. Moreover, the model parameters are estimated per day using all strikes and maturities to consider the stability of the estimates. The number of unique strikes in the data set is equal to sixteen and the number of unique maturity dates equals seven. Since not every strike is issued and/or quoted for all maturities, the number of market prices per date corresponds to 74. In addition, in our estimation routine, we will use the optimal affine parameters at time t as starting values for the estimation at time t + 1.

(33)

7. Results

The estimation routine is implemented in MatLab. The estimation method is twofold. First, we estimate the model parameters on the full set of December 2015 AEX call option prices. Second, we estimate the model parameters for each December 2015 day separately. The estimation is performed for both models as discussed in section 5.5. In the first model, it is assumed that the logarithm of the stock price follows a Vasicek process. While in the second model, it is assumed that the logarithm of the stock price follows a Heston process. Hence, stochastic volatility is added. The results for both models are presented in this chapter.

7.1. Results for the full set of market prices

The estimation results related to the Vasicek model are depicted in Table 7.1. The results related to the Heston model are depicted in Table 7.2. In order to be able to compare the parameters amongst the models, the results are presented in terms of the model parameters and not in terms of the affine parameters. The tables show the optimal parameters based on all December 2015 AEX call option prices, i.e. 1628 prices. The last column presents the minimized objective value based on equation (6.2). The objective function value is shown to compare the quality of the estimation amongst the models.

Model σ β b L

Vasicek model 0.1813 0.0082 -0.0838 4,559.53

Table 7.1.: The optimal Vasicek parameters for the full December 2015 set. The optimal pa-rameters are presented in terms of the model papa-rameters, refer to section 3.2.

Model µ ρ κ θ σv L

Heston model -0.0209 -1.0000 19.2740 0.0200 0.4396 2,671.53

Table 7.2.: The optimal Heston parameters for the full December 2015 set. The optimal pa-rameters are presented in terms of the model papa-rameters, refer to section 3.3.

The results in Table 7.1 and 7.2 indicate that using the estimated parameters of the Heston process for modelling option prices result in a smaller aggregated squared differ-ence between market and model prices compared to the Vasicek model. The additional model parameters of Heston, i.e. five versus three parameters, induce more degrees of freedom. Hence, the model is better capable of describing the dynamics of the underly-ing process of the option prices. The additional stochastic volatility process in case of

(34)

Heston, enables the model to include the dynamics of the volatility itself. The estima-tion results can be used to specify a model for Monte Carlo simulaestima-tion purposes. Recall from section 5.5 that the processes are defined such that the dynamics of log S(t) are modeled, respectively d log S(t). Based on the parameter estimates of Table 7.1 and 7.2, the following Monte Carlo specifications for Vasicek and Heston result:

Vasicek model: d log S(t) = (−0.0838 − 0.0082 log S(t))dt + 0.1813dW (t)

Heston: d log S(t) = (−0.0209 − 0.5000v(t))dt −pv(t)dW1(t)

dv(t) = 19.2740(0.0200 − v(t))dt + 0.4396pv(t)dW1(t)

Note that the specifications are in terms of continuous time. For financial applications, one often discretizes the processes using ∆ log S(t) = log S(t + ∆) − log S(t). Note that ∆ represents a discrete time step. The smaller ∆, the better the approximation.

7.2. Estimation results for the Heston model specification

Although based on a relatively large number of market prices, the previous section presents a single estimate for the model parameters. In this section, the parameters for Vasicek and Heston are estimated for each December 2015 date separately. The daily estimates allows one to compare the stability of the parameters. The drawback of course, is that the estimation is based on less observations, i.e. 74 market prices. The results over the month December are presented in the Figures 7.2. Moreover, some elementary statistics are presented in Table 7.3.

Standard

Parameter Mean deviation Max Min

Vasicek: σ 0.1803 0.0057 0.1729 0.1973 Vasicek: β 0.0144 0.0173 -0.0224 0.0443 Vasicek: b -0.1219 0.1051 -0.3044 0.0981 Heston: µ -0.0211 0.0011 -0.0226 -0.0190 Heston: ρ -0.9897 0.0237 -1.0000 -0.9018 Heston: κ 18.2803 6.9277 8.2868 33.0635 Heston: θ 0.0204 0.0010 0.0189 0.0221 Heston: σv 0.4166 0.0910 0.2449 0.5928

Table 7.3.: Statistics of the daily parameter estimates for AEX December 2015 call options.

The figures and statistics show the optimal parameters based on the daily market prices. For the Vasicek model the volatility of the process σ shows stable results with a relatively small standard deviation. The mean-reversion speed β and the long-term average b are more volatile. Note that for the stable parameters, the estimate based on the full set of market prices is in accordance with the mean. Refer to Table 7.1 and 7.3.

(35)

For the Heston model, the average µ, correlation parameter ρ, θ and the volatility of the volatility parameter σv show stable results with a relatively small standard deviation.

The κ is somewhat more volatile. Note the similarity in the patterns between the κ and σv estimates. Moreover, the average correlation parameter is close to its boundary

-1.0000.

Table 7.2 already indicated that based on all market prices, the correlation parame-ter of Heston is estimated as -1.0000. In case the correlation parameparame-ter is -1.0000, the model reduces to a form that is driven by a single Brownian motion. The resulting pro-cess can be rewritten as a one-factor model that is driven by a single Brownian motion that depends on v(t). The estimate makes one wonder if a two-factor model is indeed necessary for this set of market prices. Of course, the estimation results are for December 2015 only. Further research should test the hypothesis related to the number of factors needed. In order to provide some initial empirical evidence, we have also estimated the model parameters on AEX call option prices for other months. In addition to December 2015, we have gathered data for earlier months. The results are presented in Appendix A.9. Note that the results are in accordance with the December 2015 results.

(36)

Figure 7.2.: The graphs show the parameters for Vasicek and Heston based on the daily AEX option prices for December 2015. The top left, top right and upper middle left figures correspond respectively to σ, β and b of Vasicek. The upper middle right, lower middle left, lower middle right, lower left and lower right figures correspond respectively to µ, ρ, κ, θ and σv of Heston. The red lines are the estimates for the

full December 2015 set from Table 7.1 and 7.2.

Of course, we are also interested in the actual estimated prices. Therefore, we plot seven figures that show the estimated option prices for both Vasicek and Heston against the market prices. The graphs corresponds to December 9, 2015 and the seven different maturities. The strikes are on the horizontal axis and prices on the vertical axis. The figures are added to the Appendix in Figure A.8.

From the results, we conclude that the Heston model provides a better fit for the call option prices compared to Vasicek. The two additional parameters enable the Heston model to better capture the underlying dynamics. We can also observe that for op-tions with a larger expiry and opop-tions further away from the strike, a relatively larger performance difference between Heston and Vasicek is observed.