Hlavní práce71408_tomp03.pdf, 1.3 MB Stáhnout

(1)

Fakulta financí a účetnictví katedra bankovnictví a pojišťovnictví

studijní obor: Finanční inženýrství

Vliv intradenní volatility na ocenění opcí

Autor bakalářské práce: Bc. Pavel Tomek

Vedoucí bakalářské práce: Ing. Milan Fičura, Ph.D.

Rok obhajoby: 2020

(2)

Čestné prohlášení

Prohlašuji, že bakalářskou práci na téma „Vliv intradenní volatity na ocenění opcí“ jsem vypracoval samostatně a veškerou použitou literaturu a další prameny jsem řádně označil a uvedl v přiloženém seznamu.

V Praze dne 21. 8. 2020

(3)

Acknowledgement

I would like to hereby express my gratitude to Ing. Milan Fičura, Ph.D. for his assistance and guidance during writing of this thesis.

(4)

This thesis provides an overview of existing theory of both introductory and more advanced volatility modelling and option pricing topics. GARCH model with various extensions and stochastic volatility models under Bayesian inference are dealt with in greater detail, spline extension for intraday volatility periodicity modelling in GARCH framework is suggested. The theory is exploited to estimate eleven different volatility model specifications for Euro Stoxx 50 index, two best performing models are selected. These models are used in Monte Carlo simulations for valuation of options on Euro Stoxx 50 index under the assumption of stochastic volatility. Impact of intraday volatility on theoretical option values is studied and isolated.

Keywords: Option pricing, Intraday volatility, GARCH, Stochastic volatility models, Bayesian inference, Monte Carlo simulation.

JEL Classification: C11, C15, C32, C53, G12

Abstrakt

Tato diplomová práce podává přehled jak základní, tak pokročilejší existující teorie týkající se modelování volatility a oceňování opcí. Podrobněji je popsán model GARCH spolu s různými rozšířeními a modely stochastické volatility odhadované Bayesovskými metodami, je navrženo rozšíření GARCH modelu o spline funkci pro lepší zachycení intradenní sezónnosti volatility.

S využitím popsané teorie je odhadnuto celkem jedenáct modelů volatility s různou specifikací pro index Euro Stoxx 50, z těchto jsou vybrány dva nejlepší modely. Tyto dva modely jsou poté použity pro Monte Carlo simulace za účelem ocenění opcí na index Euro Stoxx 50 za předpokladu nekonstantní volatility. Je zkoumán a izolován vliv intradenní volatility na teoretické ocenění opcí.

Klíčová slova: Oceňování opcí, Intradenní volatilita, GARCH, Modely stochastické volatility, Bayesovské odhady, Monte Carlo simulace.

JEL klasifikace: C11, C15, C32, C53, G12

(5)

APF ... auxiliary particle filter

EWMA ... exponentially weighted moving average JPR ... Jacquier, Polson and Rossi (1994) LL ... logarithmic likelihood

LR ... likelihood ratio

LWMA ... linearly weighted moving average MC ... Monte Carlo

MCMC ... Markov chain Monte Carlo SV ... stochastic volatility

SVB ... stochastic volatility model – basic version SVC ... stochastic volatility model with correlated errors SVCT ... stochastic volatility model with correlated errors and

tails

SVCJ ... stochastic volatility model with correlated errors and jumps

SVCTJ ... stochastic volatility model with correlated errors, tails, and jumps

SX5E ... Euro Stoxx 50 index

(6)

1

Introduction

Volatility is essential to various applications in finance and financial engineering ranging from asset pricing to risk management. When looking at majority of plots that show evolution of a price of both financial or non-financial asset, one can easily observe that the plotted curve is not smooth, and its shape could be probably best described as “saw-looking”. In the most intuitive way, volatility can be thought of as a level of price variability and in turn asset return variability, which gives the curve its characteristic shape. However, volatility cannot be observed directly (we say it is a latent variable) and as such must be estimated from historical time series of asset returns, which represent one realization of the otherwise unknown data generating process.

Naïve approach would dictate to model the volatility as sample standard deviation of historical returns time series. Such approach is insufficient for vast majority of applications though because volatility is anything but constant over time. Since this is a well-known fact for a long time, many approaches to volatility modelling were developed as briefly discussed for example in Andersen and Shephard (2008).

In this thesis, we aim to provide an overview of some existing theory on (intraday) volatility modelling and option pricing. Then, our goal is to estimate models that would sufficiently capture intraday volatility process of Euro Stoxx 50 index and to apply these models to price options on the index under stochastic volatility assumption. We want to provide an evidence as to whether intraday volatility changes should theoretically impact option valuations and isolate the volatility effect from other option pricing inputs.

This thesis first briefly describes volatility in the context of financial markets and stochastic processes and discusses properties of volatility in the first chapter. The second chapter first deals with basic approaches to volatility modelling, then autoregressive conditional heteroskedasticity (“ARCH”) and stochastic volatility (in the narrow sense) frameworks are described in detail to provide the reader with a toolkit to model nearly all the volatility properties present in the hidden process.

Beginning of the third chapter provides short vocabulary concerning options and describes non-arbitrage conditions on option markets. Then option pricing basics under Black and Scholes (1973) and Merton (1973) framework are outlined and the impact of changing volatility on

(8)

3

option prices in the context of so-called Greeks is described. Lastly, we discuss the possibilities of option pricing under the relaxed assumption of constant volatility throughout the life of an option.

The last, fourth, chapter builds on the theory from the second and the third chapter. We estimate six different ARCH-type models and five stochastic volatility models that account for intraday volatility properties on various levels of complexity. For each we conduct statistical inference to select the best-performing model for each of the volatility model families. Then we apply these two models in the last part of the chapter in which we price July 2020 call options on Euro Stoxx 50 index with the use of Monte Carlo simulations.

(9)

4

1. Volatility in financial markets

1.1. Volatility and underlying stochastic processes

Trading in active financial markets could be deemed more or less continuous, which is especially true for FX markets where trading takes place 24 hours each day except for weekends and bank holidays, or for spot power markets that are open 24/7. On the other hand, historical data are usually available only for limited number of discrete time points or time intervals.

Therefore, depending on circumstances, it can be convenient to model underlying data generating process as both discrete and continuous. Discrete data generating process of asset returns can be specified as

𝑟_𝑡 = 𝑝_𝑡− 𝑝_𝑡−1 = 𝜇_{𝑡|𝑡−1}+ 𝜎_{𝑡|𝑡−1}𝜖_𝑡 (1)

where 𝑟_𝑡 represents the logarithmic return, 𝑝_𝑡 is the log-price of underlying asset, 𝜇_{𝑡|𝑡−1} is the conditional expected log-return defined as 𝑢_{𝑡|𝑡−1}= 𝐸[𝑟_𝑡|ℱ_𝑡−1] where ℱ_𝑡−1is the filtration of 𝜎-algebra that is deterministic at time 𝑡 − 1 (i.e. the realization of historical prices is known), 𝜎_{𝑡|𝑡−1} is the conditional standard deviation (conditional volatility, or conditional heteroskedasticity when we take its square) given by 𝜎_{𝑡|𝑡−1}² = 𝐷[𝑟_𝑡|ℱ_𝑡−1] = 𝐸 [(𝑟_𝑡− 𝜇_{𝑡|𝑡−1})²|ℱ_𝑡−1], and 𝜖_𝑡 ~ 𝒩(0,1) is an independent and identically distributed (“i.i.d.”) white noise where 𝒩(0,1) represents the standard normal distribution with zero mean and unit variance.

Brownian motion can be considered one of the cornerstones of financial mathematics and financial modeling, as it is probably the most used dynamics for continuous modelling of asset price behavior. It can be specified as

𝑑𝑃_𝑡 = 𝜇_𝑡𝑑𝑡 + 𝜎_𝑡𝑑𝑊_𝑡 (2)

where 𝑑𝑃_𝑡 is the instantaneous price change defined as 𝑑𝑃_𝑡= lim

Δ→0(𝑃_𝑡− 𝑃_𝑡−Δ), 𝑑𝑡 represents infinitesimal time increment defined as 𝑑𝑡 = lim

𝑁→∞(^𝑇

𝑁) where 𝑇 denotes the entire time period for which an asset is observed, and 𝑑𝑊_𝑡 is the Wiener process increment defined as 𝑑𝑊_𝑡= 𝜖_𝑡√𝑑𝑡 where 𝜖_𝑡 is 𝒩(0,1) i.i.d. white noise as above. 𝜇_𝑡 and 𝜎_𝑡 represent the drift rate and return

(10)

5

volatility. The Wiener process does not depend on its previous values (has Markov property), thus its increments are independent and

𝑊_𝑇 = ∑ 𝑑𝑊_𝑡

𝑁

𝑡=0

= ∫ 𝑑𝑊_𝑡

𝑇

0

(3)

as 𝑁 = ^𝑇

𝑑𝑡 approaches infinity. From the Equation (3) and standard normal distribution of 𝜖 follows that variance of 𝑊_𝑇 can be computed as

𝐷[𝑊_𝑇] = ∫ 𝑑𝑡

𝑇

0

= 𝑇 (4)

and standard deviation of 𝑊_𝑇 simply as √𝑇 (Hull 2012, pg. 282-284). Because the Wiener process 𝑊_𝑡 is the only source of uncertainty in Equation (2), it follows that also volatility of any asset following the Brownian motion process grows with square root of time.

However, one can observe that price volatility tends to grow in absolute terms as market price of the asset grows, which is not well captured by the Brownian motion in Equation (2). Such price dynamics are more realistically represented by the generalized Geometric Brownian motion (“GBM”) defined as

𝑑𝑃_𝑡= 𝜇_𝑡𝑃_𝑡𝑑𝑡 + 𝜎_𝑡𝑃_𝑡𝑑𝑊_𝑡 (6) or we can write

𝑑𝑃_𝑡

𝑃_𝑡 = 𝜇_𝑡𝑑𝑡 + 𝜎_𝑡𝑑𝑊_𝑡 = lim

Δ→0(𝑃_𝑡− 𝑃_𝑡−Δ

𝑃_𝑡 ) + 𝜎_𝑡𝑑𝑊_𝑡

≈ lim

Δ→0(𝑝_𝑡− 𝑝_𝑡−Δ) + 𝜎_𝑡𝑑𝑊_𝑡= 𝐸[𝑟_𝑡| ℱ_𝑡−1] + 𝜖_𝑡 (7) where 𝑝_𝑡 is the logarithm of asset price and 𝑟_𝑡 is the instantaneous log-return, 𝜖_𝑡~𝒩(0, 𝜎_𝑡²).

Equation (7) is consistent with Sharpe (1964) CAPM model which stipulates that investors require higher returns as volatility of assets (thus their riskiness) grows to be sufficiently compensated for the risks taken. Therefore, under Equations (6) and (7) market price of the assets grows more in absolute value but at more or less constant rate and similarly the price

“swings” must be more pronounced in absolute terms but more or less stable in relative terms.

Should the “absolute” volatility remain the same as the price of given asset grows, investors

(11)

6

would require lower compensation as the asset would become less risky and we would witness linear growth of the asset price instead which lead to the “basic” Brownian motion defined in Equation (2).

If 𝑝_𝑡 follows a ℱ_𝑡-adapted càdlàg¹ process with finite variation, then it can be also decomposed as

𝑝_𝑡= 𝑝₀+ 𝐴_𝑡+ 𝑀_𝑡 (8)

with 𝐴₀ = 𝑀₀ = 0, where 𝐴 is finite variance process and 𝑀 is locally square-integrable local martingale with bounded jumps (see for example Protter 2004, pg. 101-103). Andersen et al.

(2003) apply somewhat stronger assumption that 𝐴 is a finite variation and predictable mean component and further decomposition of 𝐴 and 𝑀 leads to

𝑝_𝑡 = 𝐴^𝑐_𝑡+ Δ𝐴_𝑡+ 𝑀_𝑡^𝑐+ Δ𝑀 , (9) where the predictable mean components 𝐴^𝑐 and Δ𝐴 are continuous and pure jump processes, respectively, and the local martingales 𝑀^𝑐 and Δ𝑀 are continuous sample-path and compensated jump processes, respectively. Such decomposition can be used when dealing with quadratic variation and covariation processes of returns and estimation of realized volatility.

1.2. Properties of volatility

Volatility is a latent variable, as was already mentioned in the introduction (although Andersen et al. (2003) argue otherwise to some extent) and as such needs to be estimated from historical prices and returns. Although some models assume that the level of uncertainty driving any of the underlying processes mentioned above is stable throughout the time (for example the famous Black-Scholes (1973) option pricing model), such assumption is very strong and unrealistic in most circumstances.

Volatility itself follows a stochastic process that could be, most generally, described by process of Ornstein-Uhlenbeck (“OU”) type

𝑑𝜎_𝑡² = −𝜆𝜎_𝑡²𝑑𝑡 + 𝑑𝑧_𝑡 (10)

1 From French ‘continue à droite, limite à gauche" - a process that is right continuous with left limit.

(12)

7

where 𝜆 is a (usually positive) parameter and 𝑑𝑧 is a Lévy process (that can in general be correlated with the Wiener process mentioned in Equation (2)) with positive increments. One can observe that there are time periods when “nothing happens” in the markets and the asset price evolves rather smoothly, i.e. volatility is low, and sometimes the price wildly fluctuates, i.e. volatility is high. This phenomenon is commonly termed volatility clustering and together with another important property of volatility, mean reversion, is very well captured by models of mean-reverting OU type for example described by Stein and Stein (1991) as

𝑑𝜎_𝑡² = −𝜆(𝜎_𝑡²− 𝜃)𝑑𝑡 + 𝛿𝑑𝑧_𝑡 , or more conveniently as

𝑑𝜎_𝑡² = 𝜆(𝜃 − 𝜎_𝑡²)𝑑𝑡 + 𝛿𝑑𝑧_𝑡 (11) where 𝜆, 𝜃 and 𝛿 are fixed constants, 𝜃 (a shift in the OU process) can be interpreted as the long- term level of volatility to which the volatility tends to revert – the (𝜃 − 𝜎_𝑡²) term becomes negative if volatility is higher than its mean and vice versa the term is positive if the volatility is below its mean. In general, 𝜃 could be stochastic as “normal” level of volatility in financial markets can vary through time.

Although time series of asset returns tend to display close to zero autocorrelation, volatility itself is a process with long memory (also called volatility persistence) and its autocorrelation functions tend to decline only very slowly at hyperbolic rate. The effect becomes more apparent as data sampling frequency shortens, i.e. is most pronounced in the case of high frequency intraday data.

Another property important for intraday sampling frequencies is periodicity of conditional heteroskedasticity discussed for example by Andersen and Bollerslev (1997). We can observe a U-shape of intraday volatility pattern in the case of equities for which the trading activity is high just at the beginning and at the end of trading hours. In case of currencies, the volatility pattern depends on which global financial centers are currently open for business. We can observe increased volatility from 8:00 CET when London starts trading and even higher volatility from ca. 14:00 CET when the London’s trading session overlaps with New York’s.

(13)

8

In some cases, so-called leverage effect is present in the volatility process, which means that impact of positive and negative changes of (shocks to) market prices is not symmetric (Black, 1976). For example, on stock markets exists a negative correlation between returns and volatility. In times when stock prices drop as new presumably negative information is made public, volatility increases, which could be possibly explained by the increased uncertainty regarding company’s (or entire sector/economy) future caused by the information. However, some literature suggests (see for example Bekaert and Wu (2000) for overview) that the causality can be inverse – at times of increased uncertainty expected return demanded by investors grows and so asset value has to immediately decline. In a bullish market volatility tends to decrease as we witnessed between 2010 and 2019 which was a period of unprecedently low market volatility.

For some assets, e.g. commodities, the effect can be opposite, and volatility can be positively correlated with returns. It is also consistent with growth of uncertainty because the price of commodities tends to spike in times of unexpected shortages of given commodity. Supply will likely restore at some time point in the future, thus forward prices do not have to be impacted, but spot prices grow as convenience yield increases.

Distribution of volatility of stocks is leptokurtic and positively skewed with very fat right tail for equities. By taking logarithm of the volatility we obtain approximately normal distribution (Andersen et al., 2000). Similar results can be obtained for exchange rate return volatility (Andersen et al., 2001) and can be also generalized for monthly standard deviation of returns as shown for example by French, Schwert and Stambaugh (1987).

(14)

9

2. Approaches to volatility modelling and prediction

There are two main streams of literature when it comes to approaches of volatility modelling which differ as to inputs used for their models. The first group of models relies on historical realizations of (logarithmic) returns on which econometric time series models are applied to estimate the latent volatility. Among main disadvantages of these models we could find (i) the fact, that we are given only one realization of the stochastic process, as already mentioned, thus the estimate can be biased, (ii) dependence of the outcome on correct specification of the model, and (iii) that they are sometimes computationally very demanding. One of the main advantages is the rather simple interpretation of these models. Representants of this group are for example ARCH, GARCH, ARIMA and ARFIMA building on concept of realized volatility², or stochastic volatility models. Out of these, the ARCH class of models and stochastic volatility models will be discussed in following sections in greater detail.

The other approach to volatility modelling uses prices of traded options to compute volatility.

Main advantage of these models is that they consider not only historical returns but are also forward-looking since expectations of market participants are accounted for in the option prices.

On the other hand, these models are oftentimes burdened by unrealistic assumptions such as constant volatility in the case of implied volatility computed from option prices using Black- Scholes option pricing model, or unavailability of all option quotations necessary to calculate model-free volatility. Moreover, these models oftentimes lead to volatility overestimation because of volatility risk premium.

2.1. Moving averages

Probably the most basic concept used for volatility modelling is the use of standard deviation.

Unbiased estimator of standard deviation over 𝑁 days on day 𝑡 can be defined as

𝜎̂_𝑡² = 1

𝑁 − 1∑(𝑟_𝑡−𝑖− 𝑟)²

𝑁

𝑖=1

(12)

2 As these concepts are not dealt with in any detail in this text, see for example Andersen et al. (1999) or Andersen et al. (2003).

(15)

10

where 𝜎̂_𝑡² is the estimate of variance rate of returns and its square root, 𝜎̂_𝑡, is the estimate of standard deviation (volatility) and 𝑟 is the mean value of the returns 𝑟. Because expected change in the price of asset 𝑟 is rather small in comparison with change in standard deviation, we can abstract from it (Hull 2012, pg. 498-499). Similarly, for sufficiently large 𝑁 we can replace the denominator 𝑁 − 1 just with 𝑁 and the estimator will be nearly unbiased and Equation (12) reduces to

𝜎̂_𝑡² = 1

𝑁∑ 𝑟_𝑡−𝑖²

𝑁

𝑖=1

. (13)

We can then easily model volatility for each 𝑡 ∈ 𝑇 where 𝑇 is the time period of interest. Such approach is called simple moving average (“SMA”) model. One only needs to select 𝑁 appropriately to minimize the trade-off between variance of estimates and how “up-to-date” the estimated volatility will be based on desired characteristics of the model. On one hand, large 𝑁 smooths out the unwanted variance of estimates as 𝑟_𝑡² is averaged over longer time period, on the other hand too large value means that even values from distant past can influence the estimate, which can give us misleading picture about present volatility, especially if there was a large spike in returns in the past. Therefore, some values of 𝑟² influence the outcome too much even though they are not that relevant anymore.

Possible solution to this problem is to use some weighting scheme, which leads to weighted moving average (“WMA”) model. Equation (13) then can be rewritten as

𝜎̂_𝑡² = ∑ 𝛼_𝑖𝑟_𝑡−𝑖²

𝑁

𝑖=1

, (13)

where

∑ 𝛼_𝑖 = 1

𝑁

𝑖=1

, 𝛼_𝑖 ≥ 0 (14)

are the weights assigned to each observation of 𝑟². In theory, any weights that fulfill the condition (14) could be used. In practice, linearly decreasing weights (linear weighted moving average, “LWMA”) or exponentially decreasing weights (exponential weighted moving

(16)

11

average, “EWMA”) are used as we move back through time. For LWMA then (13) simplifies to

𝜎̂_𝑡² = 2

𝑁(𝑁 + 1)∑(𝑁 − 𝑖 + 1)𝑟_𝑡−𝑖²

𝑁

𝑖=1

, (15)

but under some circumstances this model insufficiently emphasizes the newest observations.

If we wish to put even larger weight on the more recent data compared to the older data, we can use the EWMA model which can better track changes in volatility. The model can be specified as

𝜎̂_𝑡² = (1 − 𝜆) ∑ 𝜆^𝑖−1𝑟_𝑡−𝑖²

𝑁

𝑖=1

, 𝜆 ∈ (0,1) , (16)

where 𝜆 defines the rate at which the weights decay as we move into the past and the term (1 − 𝜆) ensures that the weights sum to one. 𝜆 close to 1 lead to very slow and nearly linear decay of weights, while at for example 𝜆 = 0.9 the term is close to zero already for 𝑁 about 80.

Theoretically, 𝑁 should extend infinitely into past which, of course, is not feasible. Luckily, we can choose 𝑁 such that 𝜆^𝑁−1 becomes very small so we can neglect it and, therefore, we can also omit the term for any higher 𝑁. EWMA model is for example used by RiskMetrics methodology developed by J. P. Morgan in 1994, with 𝜆 = 0.94 used for volatility estimation (Hull 2012, pg. 500-501).

2.2. ARCH and GARCH models

One of the major breakthroughs in the area of volatility modeling and estimation was the invention of autoregressive conditional heteroskedasticity (“ARCH”) class of models by Engle (1982). The model assumes that a random variable, here the asset log returns 𝑟, follows the process

𝑟_𝑡= 𝜇_{𝑡|𝑡−1}+ 𝜖_𝑡𝜎_{𝑡|𝑡−1}, (17)

where 𝐷[𝜖] = 1, which, if we add assumption of normality, i.e. 𝜖 ~ 𝒩(0,1) and zero conditional expected logarithmic return, is exactly Equation (1). Conditional variance 𝜎_{𝑡|𝑡−1}² then follows the autoregressive process

(17)

12

𝜎_{𝑡|𝑡−1}² = 𝜔 + ∑ 𝛼_𝑖𝑟_𝑡−𝑖²

𝑞

𝑖=1

, (18)

where 𝜔 can be interpreted as long-term variance rate and 𝛼_𝑖 are autoregression coefficients up to the lag 𝑞. Such model is called ARCH (q). We can then express 𝑟_𝑡 in terms of ℱ_𝑡−1, the information set available at 𝑡 − 1, as

𝑟_𝑡|ℱ_𝑡−1~𝒩(0, 𝜎_{𝑡|𝑡−1}² ) (19) and the variance function as

𝜎_{𝑡|𝑡−1}² = ℎ(𝑟_𝑡−1, 𝑟_𝑡−2, … , 𝑟_𝑡−𝑞, 𝜽) , 𝜽 = 𝛼₁, 𝛼₂, … , 𝛼_𝑞, 𝜔. (20) If we assume that the conditional expected value of 𝑟_𝑡 is non-zero, it can be expressed as linear combination of exogenous and lagged endogenous variables 𝒙_𝒕 from information set ℱ_𝑡−1 with regression coefficients 𝒃. We can then write

𝑟_𝑡|ℱ_𝑡−1~ 𝑁(𝒙_𝒕𝒃, 𝜎_{𝑡|𝑡−1}² ) , (21) 𝜎_{𝑡|𝑡−1}² = ℎ(𝜖_𝑡−1² , 𝜖_𝑡−2² , … , 𝜖_𝑡−𝑞² , 𝜽) , (22) where 𝜖_𝑡 = 𝑟_𝑡− 𝒙_𝒕𝒃̂ and also now 𝜖_𝑡|ℱ_𝑡−1~𝑁(0, 𝜎_{𝑡|𝑡−1}² ). We can then initially estimate 𝒃̂ for example by maximum likelihood method (typically by model of ARMA type) and obtain the residuals 𝜖. From the estimate of 𝜖 an efficient estimate of 𝜽̂ can be obtained by maximalization of log-likelihood function set as

𝑙 = 1 𝑇∑ 𝑙_𝑡

𝑇

𝑡=1

, (23)

where 𝑙_𝑡 are log-likelihood functions estimates for each observation 𝑡 ∈ 𝑇 𝑙_𝑡 = −1

2𝜎_{𝑡|𝑡−1}² − 𝜖_𝑡²

2𝜎_{𝑡|𝑡−1}² (24)

for which 𝜖, 𝒃̂ and 𝜽̂ are computed iteratively by switching between calculations of 𝜖, 𝒃̂ and 𝜽̂. If the process for 𝑟 can be assumed to have zero mean, we can compute the arch process directly

(18)

13

without the iterative switching. The ARCH(q) process with 𝜔 > 0 and 𝛼_𝑖, … , 𝛼_𝑞≥ 0 is covariance stationary, if ∑^𝑝_𝑖=1𝛼_𝑖 < 1.

The advantage of ARCH model is that it can quite well approximate reality as it recognizes that conditional variance varies through time due to the autoregressive component. At the same time, the model has in-built mean reversion feature and the variance decays to level of long-term (unconditional) variance 𝜔 when 𝜖 assume rather small values.

However, ARCH models sometimes are insufficient as to its ability to capture long memory property of conditional volatility, or too many lags are needed to accurately account for the entire autocorrelation structure. Another improvement for volatility estimation and prediction is represented by generalized autoregressive conditional heteroskedasticity (“GARCH”) model introduced by Bollerlsev (1986), which can be thought of as an ARCH model with AR (p) component extension. Then GARCH (p, q) can be written as

𝜎_{𝑡|𝑡−1}² = 𝜔 + ∑ 𝛼_𝑖𝜖_𝑡−𝑖²

𝑞

𝑖=1

+ ∑ 𝛽_𝑖𝜎_{𝑡−𝑖|𝑡−𝑖−1}²

𝑝

𝑖=1

(25)

where 𝜖_𝑡 = 𝑟_𝑡− 𝒙_𝒕𝒃̂ and as in case of Equation (22) 𝜖_𝑡|ℱ_𝑡−1~𝒩(0, 𝜎_{𝑡|𝑡−1}² ). The model predicts current conditional volatility based on unconditional long-term volatility, previous conditional volatility, and historical squared residuals. GARCH model can be similarly to ARCH estimated by maximum likelihood estimation by iterative switching between calculation of 𝜖, 𝒃̂, 𝜽̂ and 𝝈̂^𝟐, where 𝜽 = (𝜔, 𝛼₁, 𝛼₂, … , 𝛼_𝑞, 𝛽₁, 𝛽₂, … , 𝛽_𝑝). The model is stationary if and only if ∑^𝑞_𝑖=1𝛼_𝑖+ ∑^𝑝_𝑖=1𝛽_𝑖 < 1.

GARCH (1, 1) ranks among the most popular models for volatility prediction. It can be also specified as

𝜎_{𝑡|𝑡−1}² = 𝛾𝜎²+ 𝛼𝜖_𝑡−1+ 𝛽𝜎_{𝑡−1|𝑡−1} (26)

where 𝜎² is the long-term variance rate and 𝛾 its weight. Also 𝛾 + 𝛼 + 𝛽 = 1 since the weights must sum to unity. Even though the model directly depends only on 𝜖_𝑡−1, Bollerslev (1986) shows that the lag distribution decays at quite slow exponential rate due to the AR (1) term

(19)

14

which indicates that the model deals with long memory quite well (yet it can be improved as shown later).

When we set 𝛾 = 0, 𝛼 = (1 − 𝜆) and 𝛽 = 𝜆, equation (26) simplifies to EWMA model

𝜎_{𝑡|𝑡−1}² = (1 − 𝜆)𝜖_𝑡−1² + 𝜆𝜎_{𝑡−1|𝑡−1}² . (27) It can be shown that the last term 𝜆𝜎_{𝑡−1|𝑡−1}² converges to zero as we recursively substitute for 𝜎². At some point it is reasonable to omit the term (based on rate of decline of 𝜆) and Equation (27) reduces just to Equation (16).

Similarly, it can be shown that GARCH (p, q) can be transformed to infinite dimensional ARCH (∞) process as

𝜎_{𝑡|𝑡−1}² = 𝜔 (1 − ∑ 𝛽_𝑖

𝑝

𝑖=1

) + ∑ 𝛿_𝑖𝜖_𝑡−𝑖²

∞

𝑖=1

(28)

where 𝛿_𝑖 = 𝛼_𝑖 + ∑^𝑛_𝑗=1𝛽_𝑗𝛿_𝑗−1, 𝛼_𝑖 = 0 for 𝑖 > 𝑞 and 𝑛 = min(𝑝, 𝑖 − 1). Therefore, GARCH (p, q) can be approximated by ARCH (Q), where Q is set arbitrarily to obtain desired level of accuracy.

Although the GARCH model approximates reality reasonably well, various GARCH-based alternatives were developed to further enhance properties of the model. One such alternative is integrated GARCH (“IGARCH”) model introduced by Engle and Bollerslev (1986) to better account for long memory. The model is of the same specification as GARCH (p, q) in Equation (25) and is integrated of order 𝑑 = 1 if ∑^𝑞_𝑖=1𝛼_𝑖 + ∑^𝑝_𝑖=1𝛽_𝑖 = 1. Therefore, shocks presented to the system today are relevant for all future predictions of conditional variance, i.e. variance becomes persistent.

If we use a lag operator and denote 𝐴(𝐿) = (𝛼₁𝐿 + 𝛼₂𝐿² + ⋯ + 𝛼_𝑞𝐿^𝑞) and 𝐵(𝐿) = (𝛽₁𝐿 + 𝛽₂𝐿²+ ⋯ + 𝛽_𝑝𝐿^𝑝) where 𝐿𝑥_𝑡= 𝑥_𝑡−1, we can alternatively specify GARCH (p, q) as

𝜎_{𝑡|𝑡−1}² = 𝜔 + 𝐴(𝐿)𝜖_𝑡−1² + 𝐵(𝐿)𝜎_{𝑡|𝑡−1}² . (29) Then we can rearrange the equation to get ARMA (m, p) process of 𝜖²,

(1 − 𝐴(𝐿) − 𝐵(𝐿))𝜖_𝑡² = 𝜔 + (1 − 𝐵(𝐿))𝑣_𝑡 (30)

(20)

15

where 𝑣_𝑡 ≡ 𝜖_𝑡²− 𝜎_{𝑡|𝑡−1}² . When (1 − 𝐴(𝐿) − 𝐵(𝐿)) contains unit root we get IGARCH (p, q) model mentioned above.

If we denote 𝜙(𝐿) = (1 − 𝐴(𝐿) − 𝐵(𝐿))(1 − 𝐿)⁻¹ we can rewrite Equation (30) also as 𝜙(𝐿)(1 − 𝐿)^𝑑𝜖_𝑡² = 𝜔 + (1 − 𝐵(𝐿))𝑣_𝑡 (31) which corresponds with IGARCH (p, q) model if 𝑑 = 1. For 0 < 𝑑 ≤ 1 we get a fractionally integrated³ GARCH, FIGARCH (p, d, q), model described by Baillie, Bollerslev and Mikkelsen (1986). (1 − 𝐿)^𝑑 binomial expansion can be conveniently expressed in terms of hypergeometric function

(1 − 𝐿)^𝑑 = 𝐹(−𝑑, 1, 1; 𝐿) = Γ(−𝑑)⁻¹∑ Γ(𝑘 − 𝑑)Γ(𝑘 + 1)⁻¹𝐿^𝑘

∞

𝑘=0

, (32)

where 𝐹(⋅) is the hypergeometric function⁴ and Γ(⋅) represents the gamma function. To directly express function for conditional variance 𝜎_{𝑡|𝑡−1}² , we can further rewrite Equation (31) as

𝜎_{𝑡|𝑡−1}² = 𝜔(1 − 𝐵(𝐿))⁻¹+ [1 − (1 − 𝐵(𝐿))𝜙(𝐿)(1 − 𝐿)^𝑑]𝜖_𝑡²

≡ 𝜔(1 − 𝐵(𝐿))⁻¹+ 𝜆(𝐿)𝜖_𝑡² , (33) where 𝜆(𝐿) = 𝜆₁𝐿 + 𝜆₂𝐿²+ ⋯. Therefore, it can be seen, that ARCH representation can be yet again used to denote FIGARCH (p, d, q) with all 𝜆_𝑘≥ 0 for 𝑘 = 1, 2, … in order to assure that conditional variance is positive almost surely for all 𝑡.

As already mentioned, for 𝑑 = 1 FIGARCH equals IGARCH, for which unit shocks 𝛾(1) introduced in the system converge to 𝛾(1) = (1 + 𝐿)(1 − 𝐵(1))(1 − 𝐴(1) − 𝐵(1))⁻¹ where 𝐴(1) = ∑^𝑞_𝑖=1𝛼_𝑖 and 𝐵(1) = ∑^𝑝_𝑖=1𝛽_𝑖. For FIGARCH with 0 < 𝑑 ≤ 1 the hyperbolic decay rate will eventually dominate the long memory process dictated by 𝜆_𝑘 thus 𝛾(1) will slowly converge to 0. On the other hand, for 𝑑 > 1 we will get an explosive process and the conditional

3 See for example Arlt and Arltova (2009) for more details about fractional integration.

4 The hypergeometric function is formally defined as

𝐹(𝑚, 𝑛, 𝑠; 𝑥) = Γ(𝑠)Γ(𝑚)⁻¹Γ(𝑛)⁻¹∑^∞_𝑘=0Γ(𝑚 + 𝑘)Γ(𝑛 + 𝑘)Γ(𝑠 + 𝑘)⁻¹Γ(𝑘 + 1)⁻¹𝑥^𝑘, which for 𝑚 = −𝑑 and 𝑥 = 𝐿 reduces to the equation (32).

(21)

16

variance will grow to infinity. Such behavior of variance is not reasonable though, and, therefore, not too useful in practice. Setting 𝑑 = 0, we get the regular GARCH process for which 𝛾(1) dies out at faster exponential rate than in the case of FIGARCH. Unlike GARCH, FIGARCH is not weakly stationary, but it can be shown that it is ergodic and strictly stationary for 0 < 𝑑 ≤ 1 (Baillie, Bollerslev and Mikkelsen 1986).

Previously we discussed that volatility reacts asymmetrically to positive and negative shocks.

One of models that deals with this property is exponential GARCH (“EGARCH”) introduced by Nelson (1991) defined as

ln(𝜎_{𝑡|𝑡−1}² ) = 𝜔_𝑡+ ∑ 𝛼_𝑖𝑔(𝜖_𝑡−𝑖)

∞

𝑖=1

, 𝛼₁ = 1 (34)

or

ln(𝜎_{𝑡|𝑡−1}² ) = 𝜔_𝑡+ ∑ 𝛼_𝑖𝑔(𝜖_𝑡−𝑖) + ∑ 𝛽_𝑖ln(𝜎_{𝑡−1|𝑡−𝑖−1}² )

𝑝

𝑖 𝑞

𝑖

, 𝛼₁ = 1 (35)

where 𝜖 is i.i.d. with mean zero and unit variance. 𝜔_𝑡= 𝜔 + ln(1 + 𝑁_𝑡𝛿) is interpreted as unconditional variance for day 𝑡 with the second term representing adjustment for non-trading days. 0 < 𝛿 ≤ 1 is an adjustment parameter and 𝑁_𝑡 is number of non-trading days between 𝑡 − 1 and 𝑡. Furthermore

𝑔(𝜖_𝑡) ≡ 𝜃𝜖_𝑡+ 𝛾(|𝜖_𝑡| − 𝐸|𝜖_𝑡|), (36) which allows the conditional variance to depend asymmetrically on positive and negative shocks. Assuming that 𝛾 > 0 and 𝜃 < 0, i.e. the leverage effect is present, the curve is steeper over the interval −∞ < 𝜖_𝑡 ≤ 0 where the both terms take the same sign and “flatter” over the interval 0 < 𝜖_𝑡 < ∞ for which the terms assume opposite signs and hence at least partially offset. The parametrization in terms of logarithmic variance brings the advantage that the log-process can assume any value yet the conditional variance itself remains positive. Therefore, non-negativity constraint on model parameters does not apply.

Another possibility to add leverage effect into GARCH is Engle’s (1990) model called asymmetric GARCH (“AGARCH”) which can be specified as

(22)

17

𝜎_{𝑡|𝑡−1}² = 𝜔 + 𝛼|𝜖_𝑡−1|^𝑑 + 𝛾𝜖_𝑡−1+ 𝛽𝜎_{𝑡−1|𝑡−2}² (37) where 𝑑 is an unknown but constant power of absolute lagged residuals, 𝜔, 𝛼, 𝛾 and 𝛽 are the model parameters. If we complete the square by assuming 𝑑 = 2 the model can be also rewritten as

𝜎_{𝑡|𝑡−1}² = 𝜔^′+ 𝛼(𝜖_𝑡−1+ 𝛾^′)² + 𝛽𝜎_{𝑡−1|𝑡−2}² (38) where 𝜔^′= 𝜔 −^𝛾²

4𝛼 and 𝛾^′= ^𝛾

2𝛼. From representation of AGARCH from Equation (37) follows that for 𝛾 < 0 the impact of positive shocks on volatility is reduced as the 𝛼 and 𝛾 terms take opposite signs if 𝛼 > 0, while for negative shocks the effect is magnified. For 𝛾 > 0 the opposite is true. Equation (38) shows that AGARCH with 𝑑 = 2 is in fact only shifted version of GARCH as conditional variance is minimized for 𝜖_𝑡−1= −𝛾′ instead of 0. Then the leverage effect is accounted for only seemingly as based on sign of 𝛾 the true result of AGARCH specification leads only to the shift to one or the other direction but not to truly different sensitivity to negative and positive shocks.

Inclusion of interaction with 𝜎_{𝑡−1|𝑡−1} into the 𝛼 term leads to nonlinear AGARCH (“NAGARCH”)

𝜎_{𝑡|𝑡−1}² = 𝜔 + 𝛽𝜎_{𝑡−1|𝑡−2}² + 𝛼(𝜖_𝑡−1+ 𝛾𝜎_{𝑡−1|𝑡−2})² (39) or VGARCH

𝜎_{𝑡|𝑡−1}² = 𝜔 + 𝛽𝜎_{𝑡−1|𝑡−2}² + 𝛼(𝜖_𝑡−1𝜎_{𝑡−1|𝑡−2}⁻¹ + 𝛾)² (40) which were introduced and are further discussed by Engle and Ng (1993).

Yet another specification that leads to introduction of leverage effect into GARCH model was developed by Glosten, Jagannathan and Runkle (1993), hence called “GJR-GARCH”. They address the effect directly by adding an interaction between indicator dummy variable 𝐼 which assumes value of 1 for 𝜖_𝑡> 0 into GARCH (1, 1) which leads to

𝜎_{𝑡|𝑡−1}² = 𝜔 + 𝛼𝜖_𝑡−1² + 𝛾𝐼(𝜖_𝑡−1> 0)𝜖_𝑡−1² + 𝛽𝜎_{𝑡−1|𝑡−2 .}² (40)

(23)

18

For comparison of some of the models with leverage we generated a time series of 1 000 returns under the assumption of random walk process with normally distributed increments 𝜖_𝑡= 𝒩(0, 𝜎_{𝑡|𝑡−1}), where 𝜎_{𝑡|𝑡−1}= exp(ℎ_{𝑡|𝑡−1}/2) for which

ℎ_𝑡 = 𝜇 + 𝛽(𝜇 − ℎ_𝑡−1) + 𝐽_𝑡𝑍_𝑡^ℎ+ 𝜈, 𝜈~𝒩(0,0.0625).

𝜇 = −12, 𝐽_𝑡 represents the jump time and assumes value of one with 1% probability and 𝑍_𝑡^ℎ~𝒩(0.1,0.04) is the jump size. Parameters of the estimated models are presented in Table 1; Figure 1 compares impact of residuals under various model specifications.

Figure 1: Comparison of GARCH, AGARCH, NAGARCH, VGARCH and GJR-GARCH Source: Computations and graphics own

AGARCH model very much differs from the other specifications as the impact of residuals is much bigger compared to the other models and very asymmetric. Also Figure 1 very well shows that GARCH and GJR-GARCH models “behave” very much the same in the case of positive residuals, but impact of negative residuals is considerably larger for GJR-GARCH.

Furthermore, we can observe the small difference between the impact of positive and negative

(24)

19

residuals in the case of the non-linear models, NAGARCH and VGARCH. Moreover, VGARCH shows strong dependency on previous volatility, hence the “wobbly” shape of the curve.

Table 1: Specification of GARCH-type models with leverage

GARCH GJR-

GARCH AGARCH NAGARCH VGARCH

𝜔 0.00000010 0.00000010 0.00000037 0.00000037 0.000000086

α 0.14209 0.13828 0.25104 0.16568 0.0000013

𝛽 0.85702 0.85638 0.81255 0.81234 0.87968

𝛾 0.03829 -0.00079 0.06051 -0.00141

𝑑 1.90672

Note: Computations own.

To deal with (intraday) volatility periodicity, Bollerslev and Ghysels (1996) developed periodic GARCH (“PGARCH”) which allows model parameters to cyclically vary through time. If we define 𝑠 as the length of one cycle in number of time units identical with sampling frequency, the process may be defined as

𝜎_{𝑡|𝑡−1}² = 𝜔_𝑠(𝑡)+ ∑ 𝛼_{𝑖𝑠(𝑡)}𝜖_𝑡−1²

𝑞

𝑖=1

+ ∑ 𝛽_{𝑖𝑠(𝑡)}𝜎_{𝑡−𝑖|𝑡−𝑖−1}²

𝑝

𝑖=1

(42)

where 𝑠(𝑡) identifies the stage of periodic cycle at time 𝑡, or specified as ARMA process

𝜖_𝑡² = 𝜔_𝑠(𝑡)+ ∑ (𝛼_{𝑖𝑠(𝑡)}+ 𝛽_{𝑖𝑠(𝑡)})

max(𝑝,𝑞)

𝑖=1

𝜖_𝑡−𝑖² − ∑ 𝛽_{𝑖𝑠(𝑡)}𝑣_𝑡−𝑖

𝑝

𝑖=1

+ 𝑣_𝑡 (43)

where 𝑣_𝑡= 𝜎_{𝑡|𝑡−1}² − 𝜖_𝑡² and 𝛼_𝑖 = 0 for 𝑖 > 𝑞 and 𝛽_𝑖 = 0 for 𝑖 > 𝑝.

The fact that the model can easily become overparametrized ranks amongst its main disadvantages. As we add more parameters PGARCH loses degrees of freedom at the rate of 𝑠(𝑝 + 𝑞 + 3). Not only for this reason, it can be reasonable in practice to impose periodic behavior only on 𝛼_{𝑖(𝑠)𝑡} coefficient, and keep all 𝛽_{𝑖𝑠(𝑡)} = 𝛽_𝑖 constant across the entire cycle as it can be realistically assumed that only market participants’ reaction to shocks differs while the speed of reversion of volatility remains constant.

(25)

20

With the use of aforementioned models, PGARCH and EGARCH, we can build fractionally integrated periodic exponential GARCH (“FI-PEGARCH”) model presented by Rossi and Fantazzini (2014). First putting together PGARCH and EGARCH we get PEGARCH for logarithmic conditional variance,

ln(𝜎_{𝑡|𝑡−1}² ) = 𝜔_𝑠(𝑡)+ ∑ 𝛼_{𝑖𝑠(𝑡)}𝑔_𝑠(𝑡)(𝜖_𝑡−𝑖) + ∑ 𝛽_{𝑖𝑠(𝑡)}ln(𝜎_{𝑡−1|𝑡−𝑖−1}² )

𝑝

𝑖=1 𝑞

𝑖=1

, 𝛼₁ = 1, (44)

where 𝑔_𝑠(𝑡)(𝜖_𝑡) ≡ 𝜃_𝑠(𝑡)𝜖_𝑡+ 𝛾_𝑠(𝑡)(|𝜖_𝑡| − 𝐸|𝜖_𝑡|) Then assuming that each period is long-range dependent we can combine Equation (45) with FIGARCH and we can further rearrange the equation as

ln(𝜎_{𝑡|𝑡−1}² ) = 𝜔_𝑠(𝑡)(1 − 𝐵_𝑠(𝑡)(1)) + 𝐵_𝑠(𝑡)(𝐿)𝜎_{𝑡−1|𝑡−2}² + (1 − 𝐿)^−𝑑(1 + 𝐴_𝑠(𝑡)(𝐿))𝑣_𝑡−1 (45) where 𝑣_𝑡 = 𝑔_𝑠(𝑡)(𝜖_𝑡) − ln(𝜎_{𝑡|𝑡−1}² ). Eventually, we can express FI-PEGARCH as infinite dimensional ARCH process

ln(𝜎_{𝑡|𝑡−1}² ) = 𝜔_𝑠(𝑡)+ (1 − 𝐿)^−𝑑(1 − 𝐵_𝑠(𝑡)(𝐿))⁻¹(1 + 𝐴_𝑠(𝑡)(𝐿)) 𝑔_𝑠(𝑡)(𝜖_𝑡−1). (46) The number of degrees of freedom declines rapidly as the number of estimated parameters for each period increases. Such problem cannot be dealt with easily under the EGARCH framework because 𝛼₁ = 1, which means that we also have to estimate periodic parameters for 𝜃_𝑠(𝑡) and 𝛾_𝑠(𝑡) from the function 𝑔_𝑠(𝑡)(𝜖_𝑡) if we want to impose periodicity on the asymmetric reactions to shocks. Therefore, we propose to model the asymmetry using the GJR-GARCH methodology, but to further reduce the number of estimated parameters we assume that the periodic changes of asymmetry are proportionate to periodic changes of 𝛼_𝑠(𝑡) parameters, which can be included by the interaction term 𝛾_𝑖𝛼_𝑠(𝑡). Then further focusing only on periodicity of reactions to changes in 𝜖_𝑡 we can write the model as

𝜎_{𝑡|𝑡−1}² = 𝜔 + (1 − 𝐿)^−𝑑{∑ 𝛼_{𝑖𝑠(𝑡)}[𝜖_𝑡−𝑖² + 𝛾_𝑖𝐼(𝜖_𝑡−1 > 0)𝑒_𝑡−𝑖² ]

𝑞

𝑖=1

} + ∑ 𝛽_𝑖𝜎_{𝑡−𝑖|𝑡−𝑖−1}²

𝑝

𝑖=1

(47)

It sounds reasonable to impose periodicity on the long-term variance rate 𝜔 as well, but again the number of estimated parameters would have to be high to achieve acceptable level of

(26)

21

smoothness of the intraday pattern. Engle and Rangel (2008) suggested to use exponential splines to deal with constant level of volatility mean reversion over the long term. Motivated by their research, we propose to incorporate the changing level of unconditional variance of intraday level to Equation (46) by cubic spline function

𝜎_{𝑡|𝑡−1}² = 𝜔 + ∑ 𝛿_𝑖𝑓_𝑖(𝑡_𝑠, 𝜏_𝑖)

𝐾

𝑖

+ (1 − 𝐿)^−𝑑{∑ 𝛼_{𝑖𝑠(𝑡)}[𝑒_𝑡−𝑖² + 𝛾_𝑖𝐼(𝜖_𝑡−1> 0)𝑒_𝑡−𝑖² ]

𝑞

𝑖=1

}

+ ∑ 𝛽_𝑖𝜎_{𝑡−𝑖|𝑡−𝑖−1}²

𝑝

𝑖=1

(48)

where 𝜔 can be interpreted as the constant unconditional variance across the entire daily cycle, 𝐾 is the number of cubic spline functions to be used, 𝑡_𝑠 tracks time that lapsed since beginning of given trading day (based on data granularity, e.g. there are 27 periods a day for 20-minute returns of EuroStoxx 50 discussed below, hence 𝑡_𝑠 = 1, … ,27), 𝛿_𝑖 are estimated spline function parameters and

𝑓(𝑡_𝑠, 𝜏_𝑖) = max(0, (𝑡_𝑠− 𝜏_𝑖)³) (49) where 𝜏_𝑖 denotes 𝑡_𝑠 (so-called “knot”) for which given spline function becomes applicable.

2.3. Stochastic volatility models

The theory underlying stochastic volatility (“SV”) models extends back at least to Clark (1973) who proposed to model speculative price series as subordinate stochastic process to deal with excess kurtosis and fat tails of returns. If we denote 𝑃_𝑡 as price at time 𝑡 driven by process set in Equation (1) with zero drift then we get random walk with process for price increments Δ𝑃_𝑡= 𝑃_𝑡− 𝑃_𝑡−1 that are i.i.d with variance 𝜎². Clark further establishes a process 𝑇(𝑡) with increments drawn from positive distribution with mean 𝛼 > 0, i.e. a process that constitutes trading clock. Then by substitution of 𝑇(𝑡) into 𝑃_𝑡 we get 𝑃_𝑇(𝑡) with increments Δ𝑃_𝑇(𝑡), a process that is subordinated to Δ𝑃_𝑡. 𝑇(𝑡) can be thought of as speed of price evolution and could depend on arrival of news that influence investors’ decision making. Therefore, more important, or more frequent news can lead to larger 𝑇(𝑡) and to larger volatility in turn. Since 𝑇(𝑡) is stochastic, it follows that volatility is also stochastic. Also, since central limit theorem is

(27)

22

applicable, variables that follow aforementioned process of Δ𝑃_𝑇(𝑡) are conditionally normally distributed.

In a broader sense, ARCH-type models can be considered of SV class, too. However, one important difference between these and “true” SV process, as we understand them in this section, is apparent. In case of ARCH processes, the latent volatility process is driven by the underlying price (return) process, which, although itself being stochastic, is already deterministic at the time when volatility is estimated as all the explanatory variables are already known. On the other hand, in case of SV models, volatility follows its own stochastic process possibly independent of the returns but influenced by past variance and also future variance.

One of the reasons why SV models gained on popularity only quite recently were the computational difficulties that econometricians faced early on as the likelihood function of these models is hard to evaluate. For example, Melino and Turnbull (1990) decided to evade this problem by applying method of moments (“MM”) to obtain their estimate. However, the moments selected for the computations are chosen somewhat arbitrarily and, more importantly, MM estimates do not reach Rao-Cramer lower bound thus lead to inefficient estimates compared to likelihood methods.

Jacquier, Polson and Rossi (1994) (“JPR”) suggested using Markov chain Monte Carlo (“MCMC”) approach based on Bayesian inference and mixture density hypothesis. In their simple model, an asset follows a process

𝑟_𝑡= 𝑒^ℎ^𝑡^/2𝜖_𝑡 (50)

h_t = 𝜇 + 𝛽ℎ_𝑡−1+ 𝜎_𝜈𝜈_𝑡

where ℎ_𝑡 = ln 𝜎_𝑡², 𝑟_𝑡 is the logarithmic return at time 𝑡, 𝜎_𝑡² is the conditional variance, 𝜇 the long-term volatility component, 𝛽 is the AR (1) coefficient and 𝜎_𝜈 is the volatility of conditional variance. Further 𝜖_𝑡, 𝜈_𝑡 ~ 𝒩(0, 1) and in this model 𝐸[𝜖_𝑡𝜈_𝑡] = 0 is assumed.

According to Bayes theorem (see for example Gelman et al 2014, pg. 7) we search for posterior distribution

𝑝(𝜃|𝑦) ∝ 𝑝(𝑦|𝜃)𝑝(𝜃) , (51)

(28)

23

where 𝑝(⋅) denotes the probability distribution, 𝜃 the parameter of interest and 𝑦 the observed data. Furthermore, 𝑝(𝜃) represents the prior distribution of 𝜃 (or unknown variable), which should summarize all information about given parameter (variable) that a researcher has prior making any inference. Based on settings of 𝑝(𝜃), one can put various emphasis on information they had beforehand. Rather flat (uninformative) prior probability distribution will let the observed data “speak for themselves” while strongly informative prior can materially skew the outcome of an analysis in some direction.

Hence, in context of volatility modelling, JPR suggest looking for posterior distribution given by

𝑝(𝒉, 𝜽|𝒓) = 𝑝(𝒓|𝒉)𝑝(𝒉|𝜽)𝑝(𝜽). (52)

where 𝜽 = (𝜇, 𝛽, 𝜎_𝜈)^′, 𝒓 = (𝑟₁, 𝑟₂, … , 𝑟_𝑇)^′ and 𝒉 = (ℎ₁, ℎ₂, … , 𝑦ℎ_𝑇)^′. In other applications one could break 𝑝(𝒉, 𝜽|𝒓) into marginal conditional distributions 𝑝(𝒉|𝜽, 𝒓) and 𝑝(𝜽|𝒉, 𝒓) and then sample iteratively sample 𝜽^𝑗 and 𝒉^𝑗 , 𝑗 = 0,1,2 … , from these distributions following a generic procedure:

1. Set initial values of 𝒉⁰ and 𝜽⁰. 𝒉⁰ can be set for example to vector of unconditional variances;

it is convenient to preset 𝜽⁰ to some values which the parameters can be expected to reasonably assume to speed up the convergence.

2. Set 𝑗 = 𝑗 + 1

3. Draw sample of 𝜽^𝑗 from 𝑝(𝜽|𝒉, 𝒓).

4. Draw samples of 𝒉^𝑗 from 𝑝(𝒉|𝜽, 𝒓).

5. Return to step 2 and repeat the sweep⁵.

One can then repeat the abovementioned algorithm 𝑀 times where 𝑀 is set to obtain desired level of accuracy of results. However, in the case of this model we can only efficiently sample from 𝑝(𝜽|𝒉, 𝒓).

Using standard linear models, we obtain samples for 𝜇, 𝛽|𝜎_𝜈, ℎ_𝑡 ~ 𝒩₂(𝜷̂, 𝜎_𝜈𝑨⁻¹), 𝜷̂ = (𝜇, 𝛽)^′, 𝑨 = 𝑿_𝒕^′𝑿_𝒕, where 𝑿_𝒕 is a matrix consisting of column of ones and a column of lagged ℎ_𝑡, and

5 One iteration of estimation of all values 𝜽 and 𝒉 is called a sweep.