# Markov Switching Model Analysis of Implied Volatility for Market Indexes with Applications to S&P 500 and DAX

Authors: Luca Di Persio; Samuele Vettori
Journal: Journal of Mathematics
 (2014)

---

## Abstract

We adopt aregime switching approach to study concrete financial time series with particular emphasis on their volatility characteristics considered in a space-time setting. In particular the volatility parameter is treated as an unobserved state variable whose value in time is given as the outcome of an unobserved, discrete-time and discrete-state, stochastic process represented by a suitable Markov chain. We will take into account two different approaches for inference on Markov switching models, namely, the classical approach based on the maximum likelihood techniques and the Bayesian inference method realized through a Gibbs sampling procedure. Then the classical approach shall be tested on data taken from the Standard & Poor’s 500 and the Deutsche Aktien Index series of returns in different time periods. Computations are given for a four-state switching model and obtained numerical results are put beside by explanatory graphs which report the outcomes obtained exploiting both smoothing and filtering algorithms used in the estimation/calibration procedures we proposed to infer on the switching model parameters.


---

## Body

## 1. Introduction

Many financial time series are characterized by abrupt changes in their behaviour, a phenomena that can be implied by a number of both endogeneous and exogeneous facts, often far from being forecasted. Examples of such changing factors can be represented, for example, by large financial crises, government policy and political instabilities, natural disasters, and speculative initiatives.Such phenomena have been frequently observed during last decade especially because of the worldwide financial crisis which originated during the first years of 2000 and is still running. Indeed such a crisis, often referred to as theGlobal Financial Crisis, has caused a big lack of liquidity for banks, both in USA, Europe, and many countries all over the world, resulting, for example, in a collapse of many financial institutions, a generalized downfall in stock markets, a decline consumer wealth, and an impressive growth of the Eurpean sovereign debt. The reasons behind these phenomena are rather complicated, particularly because of the high number of interconnected point of interests, each of which is driven by specific influences often linked between each other. Nevertheless there are mathematical techniques which can be used to point out some general characteristics able to synthesize some relevant informations and to give an effective help in forecasting future behaviour of certain macroquantities of particular interest.Although linear time series techniques, for example, the autoregressive (AR) model, the moving average (MA) model, and their combination (ARMA), have been successfully applied in a large number of financial applications, by their own nature they are unable to describe nonlinear dynamic patterns as in the case, for example, of asymmetry and volatility clustering. In order to overcome latter issue various approaches have been developed. Between them and with financial applications in mind, we recall the autoregressive conditional heteroskedasticity (ARCH) model of Engle, together with its generalised version (GARCH), and theregime switching (RS) model which involves multiple equations, each characterizing the behaviors of quantities of interest in different regimes, and a mechanism allowing the switch between them. Concerning the switching methods that can be considered in the RS framework, we would like to cite thethreshold autoregressive (TAR) model proposed by Tong in [1], in which regime switching is controlled by a fixed threshold, theautoregressive conditional root (ACR) model of Bec et al. (see [2]) where the regime switching between stationary and nonstationary state is controlled by a binary random variable, and its extension, namely, thefunctional coefficient autoregressive conditional root (FCACR) model, considered by Zhou and Chen in [3]. In particular in this work we aim at using the RS approach to model aforementioned types of unexpected changes by their dependence on an unobserved variable typically defined as theregime orstate. A customary way to formalize such an approach is given through the followingstate-space representation:
(1)
y
t
=
f
S
t
,
θ
ψ
,
t
=
1
,
…
,
n
,
n
∈
N
+
,
S
t
∈
Ω
,
where θ is the vector of problem parameters, ψ is the information set, thestate set
Ω, with |
Ω
|
=
M
∈
N
+, is the (finite) set of possible values which can be taken by the the state process S at time t, that is, S
t
∈
Ω, and f is a suitable function determining the value of the dependent variable y at any given time t
∈
1
,
…
,
n.The simplest type of structure considers two regimes; that is,Ω
=
{
1,2
} and at most one switch in the time series: in other words, the first m (unknown) observations relate regime 1, while the remaining n
-
m data concern regime 2. Such an approach can be generalized allowing the system to switch back and forth between the two regimes, with a certain probability. The latter is the case considered by Quandt in his paper published in 1962, where it is assumed to be theoretically possible for the system to switch between regimes every time that a new observation is generated. Note that previous hypothesis is not realistic in an economic context, since it contradicts the volatility clustering property, typical of financial time series.The best way to represent the volatility clusters phenomenon consists in assuming that the state variable follows a Markov chain and claiming that the probability of having a switch in the next time is much lower than the probability of remaining in the same economic regime. The Markovian switching mechanism was first considered by Goldfeld and Quandt in [4] and then extended by Hamilton to the case of structural changes in the parameters of an autoregressive process (see [5]). When the unobserved state variable that controls the switching mechanism follows a first-order Markov chain, the RS model is called Markov Switching Model (MSM). In particular the Markovian property of such a model implies that, given t
∈
2
,
…
,
n, the value S
t of the state variable depends only on S
t
-
1, a property that turns out to be useful to obtain a good representation of financial data, where abrupt changes occur occasionally.After Hamilton’s papers, Markov switching models have been widely applied, together with a number of alternative versions, to analyze different type of both economic and financial time series, for example, stock options behaviors, energy markets trends, and interest rates series. In what follows we shall often refer to the following, actually rather general, form for an MSM:(2)
y
t
=
f
(
S
t
,
θ
)
S
t
∈
{
1,2
,
…
,
M
}
with
probabilities
p
i
j
:
=
P
S
t
=
j
∣
S
t
-
1
=
i
,
i
,
j
=
1,2
,
…
,
M
,
t
=
2
,
…
,
n
,
where the terms p
i
j are nothing but the transition probabilities, from state i at time t
-
1 to state j at time t, which determine the stochastic dynamic of the process S
=
S
t
t
=
1
,
…
,
n.In Section2 we recall the classical approach to MSM following [5] and with respect to both serially uncorrelated and correlated data. Then, in Section 3, we first introduce basic facts related to the Bayesian inference; then we recall the Gibbs sampling technique and related Monte Carlo approximation method which is later used to infer on the MSM parameters. Section 4 is devoted to agoodness of fit (of obtained estimates for parameters of interest) analysis, while in Section 5 we describe our forecasting MSM-based procedure which is then used in Section 6 to analyse both theStandard & Poor’s 500 and theDeutsche Aktien indexes.


## 2. The Classical Approach

In this section we shall use theclassical approach (see, for example, [5]) to develop procedures which allow us to make inference on unobserved variables and parameters characterizing MSM. The main idea behind such an approach is splitin two steps: first we estimate the model’s unknown parameters by a maximum likelihood method; secondly we infer the unobserved switching variable values conditional on the parameter estimates. Along latter lines, we shall analyze two different MSM settings, namely, the case in which data areserially uncorrelated, and the case when they areautocorrelated.


### 2.1. Serially Uncorrelated Data

Let us suppose thatS
=
S
t
t
∈
{
1
,
…
,
T
}, T
∈
N
+, is a discrete time stochastic process represented by a first-order Markov chain taking value in some set Ω
=
{
1
,
…
,
M
}, with transition probabilities:
(3)
p
i
j
∶
=
P
(
S
t
+
1
=
j
∣
S
t
=
i
,
S
t
-
1
=
i
t
-
1
,
…
,
S
0
=
i
0
)
=
P
(
S
t
+
1
=
j
∣
S
t
=
i
)
,
∀
t
∈
{
2
,
…
,
T
}
.
Then the state-space model we want to study is given by
(4)
y
t
=
μ
S
t
+
ϵ
t
,
t
=
1,2
,
…
,
T
,
ϵ
t
t
∈
{
1
,
…
,
T
}
are
i
.
i
.
d
.
N
(
0
,
σ
S
t
2
)
,
μ
S
t
=
μ
1
S
1
,
t
+
⋯
+
μ
M
S
M
,
t
,
σ
S
t
=
σ
1
S
1
,
t
+
⋯
+
σ
M
S
M
,
t
,
S
t
∈
{
1
,
…
,
M
}
with
probabilities
p
i
j
,
where the variables S
k
,
t, (
k
,
t
)
∈
{
1
,
…
,
M
}
×
{
1
,
…
,
T
} are introduced in order to have a slightly compact equation for μ
S
t and σ
S
t: in particular, S
k
,
t
=
1 if S
t
=
k, otherwise S
k
,
t
=
0, which implies that under regime k, for k
∈
{
1
,
…
,
M
}, parameters of mean and variance are given by μ
k and σ
k.Let us underline that in (4) the y
t are our observed data, for example, historical returns of a stock or some index time series, and we suppose that they arelocally distributed as Gaussian random variable in the sense that occasionally jumps could occur for both the mean μ
S
t and the variance σ
S
t
2. In particular, we assume that σ
1
<
σ
2
<
⋯
<
σ
M and we want to estimate these M unobserved values for standard deviation, as well as the M values for the switching mean. Note that we could also take μ as a constant, obtaining the so calledswitching variance problem or σ as a constant, having aswitching mean problem. The first one is in general more interesting in the analysis of financial time series, since variance is usually interpreted as an indicator for the market’s volatility.Given the model described by (4), the conditional density of y
t given S
t is Gaussian with mean μ
S
t and standard deviation σ
S
t; namely, its related probability density function reads as follow:
(5)
f
y
t
∣
S
t
(
x
)
=
1
2
π
σ
S
t
2
e
-
(
x
-
μ
S
t
)
2
/
2
σ
S
t
2
,
x
∈
R
,
and we are left with the problem of estimating both the expectations μ
1
,
…
,
μ
M and the standard deviations σ
1
,
…
,
σ
M parameters; a task that is standard to solve by maximizing the associated log-likelihood function ln
⁡
L defined as follows:
(6)
ln
⁡
L
:
=
∑
t
=
1
T
‍
ln
⁡
(
f
y
t
∣
S
t
)
.
A different and more realistic scenario is the one characterized by unobserved values for S
t. In such a case, it could be possible to consider the MSM-inference problem as a two-step procedure consisting in(i)
estimating the parameters of the model by maximizing the log-likelihood function,
(ii)
making inferences on the state variableS
t,  t
=
1
,
…
,
T conditional on the estimates obtained at previous point.Depending on the amount of information we can use inferencing onS
t, we have(i)
filtered probabilities that refer to inferences about S
t conditional on information up to time t, namely, with respect to ψ
t,
(ii)
smoothed probabilities that refer to inferences about S
t conditional to the whole sample (history), ψ
T.In what follows we describe a step-by-step algorithm which allows us to resolve the filtering problem for a sample of serially uncorrelated data. In particular we slightly generalize the approach given in [6], assuming that the state variable S
t belongs to a 4-state space set Ω
:
=
{
1,2
,
3,4
}, at every time t
=
1
,
…
,
T. Despite 2-state (expansion/contraction) and 3-state (low/medium/high volatility regime) models being the usual choices, we decided to consider 4-state MSM in order to refine the analysis with respect to volatility levels around the mean. A finer analysis can be also performed, even if one has to take into account the related nonlinear computational growth. We first define the log-likelihood function at time t as
(7)
l
(
θ
,
ψ
t
)
=
l
(
μ
1
,
…
,
μ
M
,
σ
1
,
…
,
σ
M
,
ψ
t
)
,
where θ
:
=
(
μ
1
,
…
,
μ
M
,
σ
1
,
…
,
σ
M
) is the vector of parameters that we want to estimate. Let us note that the θ is updated at every iteration, since we maximize with respect to the function l
(
θ
,
ψ
t
) at every stage of our step-by-step procedure. In particular the calibrating procedure reads as follows.Inputs.(i)
Putl
(
θ
)
:
=
l
(
θ
,
ψ
0
)
=
0.
(ii)
Compute thetransition probabilities of the homogeneous Markov chain underlying the state variable; that is,
(8)
p
i
j
=
P
(
S
t
=
j
∣
S
t
-
1
=
i
)
,
i
,
j
∈
Ω
.
Since in the applications we can only count on return time series, we first have to calibrate with respect to the transition probabilitiesp
i
j:
(1)
choose 4 valuesσ
1
2
<
σ
2
2
<
σ
3
2
<
σ
4
2 (e.g., σ
1
=
0.1, σ
2
=
0.2, σ
3
=
0.3, and σ
4
=
0.4) and a positive, arbitrary small, constant δ
>
0, for example, δ
=
0.01,
(2)
compute for everyj
=
1,2
,
3,4 and for every t
=
1
,
…
,
T the values
(9)
b
j
≔
Φ
y
t
+
δ
-
Φ
y
t
-
δ
=
∫
y
t
-
δ
y
t
+
δ
‍
1
2
π
σ
j
2
exp
⁡
x
-
μ
j
2
2
σ
j
2
,
(3)
simulate a value in{
1,2
,
3,4
} for S
t at each time from the discrete probability vector
(10)
b
1
∑
j
=
1
4
‍
b
j
,
b
2
∑
j
=
1
4
‍
b
j
,
b
3
∑
j
=
1
4
‍
b
j
,
b
4
∑
j
=
1
4
‍
b
j
,
(4)
set the transition probabilitiesp
i
j just by counting the number of transition from state i to state j, for i
,
j
=
1
,
…
,
4, in order to obtain the followingtransition matrix:
(11)
P
*
=
p
11
p
12
p
13
p
14
p
21
p
22
p
23
p
24
p
31
p
32
p
33
p
34
p
41
p
42
p
43
p
44
.
(iii)
Compute thesteady-state probabilities
(12)
π
(
0
)
=
P
S
0
=
1
∣
ψ
0
,
…
,
P
S
0
=
4
∣
ψ
0
.
Let us note that, by definition, ifπ
(
t
) is a 4
×
1 vector of steady-state probabilities, then π
(
t
+
1
)
=
π
(
t
) for every t
=
1
,
…
,
T; moreover π
(
t
+
1
)
=
P
*
π
(
t
), and (see for example, [6, pag.  71]) we also have that π
(
t
)
=
(
A
T
A
)
-
1
A
T
0
4
1, where 0
4 is a 4
×
1 matrix of zeros and A
=
I
d
4
-
P
*
1
4, I
d
4 is the four dimensional identity matrix, while 1
4
=
(
1,1
,
1,1
), that is, the vector of steady-state probabilities is the last column of the matrix (
A
T
A
)
-
1
A
T.Next we perform the following steps, fort
=
1
,
…
,
T.Step 1.
The probability ofS
t conditional to information set at time t
-
1 is given by
(13)
P
S
t
=
j
∣
ψ
t
-
1
=
∑
i
=
1
4
‍
p
i
j
P
S
t
-
1
=
i
∣
ψ
t
-
1
,
j
=
1
,
…
,
4
.Step 2.
Compute the joint density ofy
t and S
t conditional to the information set ψ
t
-
1
(14)
f
y
t
,
S
t
=
j
∣
ψ
t
-
1
=
f
y
t
∣
S
t
=
j
,
ψ
t
-
1
×
P
S
t
=
j
∣
ψ
t
-
1
,
j
=
1
,
…
,
4
.
The marginal density of y
t is given by the sum of the joint density over all values of S
t
(15)
f
y
t
∣
ψ
t
-
1
=
∑
i
=
1
4
‍
f
y
t
∣
S
t
=
i
,
ψ
t
-
1
P
S
t
=
i
∣
ψ
t
-
1
.Step 3.
Update the log-likelihood function at timet in the following way:
(16)
l
θ
,
ψ
t
=
l
θ
,
ψ
t
-
1
+
ln
⁡
f
y
t
∣
ψ
t
-
1
,
and maximize l
(
θ
,
ψ
t
) with respect to θ
=
(
μ
1
,
…
,
μ
4
,
σ
1
,
…
,
σ
4
), under the condition σ
1
<
σ
2
<
σ
3
<
σ
4, to find the maximum likelihood estimator θ
^ for the next time period.Step 4.
Oncey
t is observed at the end of the tth iteration, we canupdate the probability term:
(17)
P
S
t
=
j
∣
ψ
t
=
f
y
t
S
t
=
j
,
ψ
t
-
1
P
S
t
=
j
ψ
t
-
1
∑
i
=
1
4
‍
f
y
t
∣
S
t
=
i
,
ψ
t
-
1
P
S
t
=
i
∣
ψ
t
-
1
,
hhhhhhhhhhhhhhhhhhhhhhih
j
=
1
,
…
,
4
,
where both ψ
t
=
{
ψ
t
-
1
,
y
t
} and f
(
y
t
∣
S
t
=
i
,
ψ
t
-
1
) are computed with respect to the estimator θ
^
=
(
μ
^
1
,
…
,
μ
^
4
,
σ
^
1
,
…
,
σ
^
4
).


### 2.2. Serially Correlated Data

In some cases it is possible to argue and mathematically test by, for example, the Durbin-Watson statistics or Breusch-Godfrey test, for the presence of a serial correlation (orautocorrelation) between data belonging to a certain time series of interest. Such a characteristic is often analyzed in signal processing scenario, but examples can be also found in economic, meteorological, or sociological data sets, especially in connection with autocorrelation of errors in related forecasting procedure. In particular if we suppose that the observed variable y
t linearly depends on its previous value, then we obtain afirst-order autoregressive pattern and the following state-space model applies:
(18)
y
t
-
μ
S
t
=
ϕ
(
y
t
-
1
-
μ
S
t
-
1
)
+
ϵ
t
,
t
=
1,2
,
…
,
T
,
ϵ
t
t
∈
{
1
,
…
,
T
}
are
i
.
i
.
d
.
N
(
0
,
σ
S
t
2
)
,
μ
S
t
=
μ
1
S
1
,
t
+
⋯
+
μ
M
S
M
,
t
,
σ
S
t
=
σ
1
S
1
,
t
+
⋯
+
σ
M
S
M
,
t
,
S
t
∈
1
,
…
,
M
with
probabilities
p
i
j
,
where p
i
j
:
=
P
(
S
t
=
j
∣
S
t
-
1
=
i
), i
,
j
=
1
,
…
,
M and S
k
,
t, and (
k
,
t
)
∈
{
1
,
…
,
M
}
×
{
1
,
…
,
T
} are the same variables introduced in the previous section; that is, S
k
,
t
=
1 if S
t
=
k, otherwise S
k
,
t
=
0.In this situation, if the stateS
t is known for every t
=
1
,
…
,
T, we need S
t and S
t
-
1 to compute the density of y
t conditional to past information ψ
t
-
1, indeed we have
(19)
ln
⁡
L
:
=
∑
t
=
1
T
‍
ln
⁡
f
y
t
∣
ψ
t
-
1
,
S
t
,
S
t
-
1
,
where
(20)
f
y
t
∣
ψ
t
-
1
,
S
t
,
S
t
-
1
(
x
)
=
1
2
π
σ
S
t
2
e
-
x
-
μ
S
t
-
ϕ
y
t
-
1
-
μ
S
t
-
1
2
/
2
σ
S
t
2
,
x
∈
R
.IfS
t are unobserved (and as before we assume that the state variable can take the four values {
1,2
,
3,4
}), we apply the following algorithm in order to resolve the filtering problem for a sample of serially correlated data.Inputs.(i)
Putl
(
θ
)
:
=
l
(
θ
,
ψ
0
)
=
0.
(ii)
Compute the transition probabilities(21)
p
i
j
=
P
S
t
=
j
∣
S
t
-
1
=
i
,
i
,
j
=
1,2
,
3,4
.
We apply the same trick as before, but firstly we have to estimate the parameterϕ: in order to obtain this value we can use the least square methods (see for example, [7]) that is,
(22)
ϕ
^
:
=
∑
t
=
1
T
‍
y
t
y
t
+
1
∑
t
=
1
T
‍
y
t
2
.
Then we computez
t
:
=
y
t
-
ϕ
y
t
-
1 for every t
=
1
,
…
,
T and consider the values b
j
:
=
Φ
(
z
t
+
δ
)
-
Φ
(
z
t
-
δ
) (we apply the Normal distribution function to z
t
+
δ instead of y
t
+
δ, as done before).
(iii)
Compute the steady-state probabilities(23)
π
(
0
)
=
P
S
0
=
1
∣
ψ
0
,
…
,
P
S
0
=
4
∣
ψ
0
,
taking the last column of the matrix(
A
T
A
)
-
1
A
T (see procedure in Section 2.1 for details).Next perform the following steps for t
=
1
,
…
,
T.Step 1.
Compute the probabilities ofS
t conditional to information set at time t
-
1, for j
=
1
,
…
,
4
(24)
P
S
t
=
j
∣
ψ
t
-
1
=
∑
i
=
1
4
‍
P
S
t
=
j
,
S
t
-
1
=
i
∣
ψ
t
-
1
=
∑
i
=
1
4
‍
p
i
j
P
S
t
-
1
=
i
∣
ψ
t
-
1
.Step 2.
Compute the joint density ofy
t, S
t, and S
t
-
1 given ψ
t
-
1:
(25)
f
y
t
,
S
t
,
S
t
-
1
∣
ψ
t
-
1
=
f
y
t
∣
S
t
=
j
,
S
t
-
1
=
i
,
ψ
t
-
1
×
P
S
t
=
j
,
S
t
-
1
=
i
∣
ψ
t
-
1
,
where f
(
y
t
∣
S
t
=
j
,
S
t
-
1
=
i
,
ψ
t
-
1
) is given by (20) and P
(
S
t
=
j
,
S
t
-
1
=
i
∣
ψ
t
-
1
) is computed in Step 1. The marginal density of y
t conditional on ψ
t
-
1 is obtained by summing over all values of S
t and S
t
-
1:
(26)
f
y
t
∣
ψ
t
-
1
=
∑
j
=
1
4
‍
∑
i
=
1
4
‍
f
y
t
∣
S
t
=
j
,
S
t
-
1
=
i
,
ψ
t
-
1
hhhhl
×
P
S
t
=
j
,
S
t
-
1
=
i
∣
ψ
t
-
1
.Step 3.
The log-likelihood function at timet is again
(27)
l
θ
,
ψ
t
=
l
θ
,
ψ
t
-
1
+
ln
⁡
f
y
t
∣
ψ
t
-
1
,
and it can be maximized with respect to θ
=
(
μ
1
,
…
,
μ
4
,
σ
1
,
…
,
σ
4
), under condition σ
1
<
σ
2
<
σ
3
<
σ
4, giving the maximum likelihood estimator θ
^ for the next time period.Step 4.
Update the joint probabilities ofS
t and S
t
-
1 conditional to the new information set ψ
t, using the estimator θ
^ computed in Step 3 by maximizing the log-likelihood function l
(
θ
,
ψ
t
):
(28)
P
S
t
=
j
,
S
t
-
1
=
i
∣
ψ
t
=
f
y
t
,
S
t
=
j
,
S
t
-
1
=
i
∣
ψ
t
-
1
f
y
t
∣
ψ
t
-
1
,
i
,
j
=
1
,
…
,
4
.
Then compute the updated probabilities ofS
t given ψ
t by summing the joint probabilities over S
t
-
1 as follows:
(29)
P
S
t
=
j
∣
ψ
t
=
∑
i
=
1
4
‍
P
S
t
=
j
,
S
t
-
1
=
i
∣
ψ
t
,
∀
j
=
1
,
…
,
4
.The Smoothing Algorithm. Once we have run this procedure, we are provided with thefiltered probabilities, that is, the values P
(
S
t
=
j
∣
ψ
t
) for j
=
1
,
…
,
4 and for each t
=
1
,
…
,
T (in addition to the estimator θ
^).Sometimes it is required to estimate probabilities ofS
t given the whole sample information; that is,
(30)
P
S
t
=
j
∣
ψ
T
=
P
S
t
=
j
∣
y
1
,
…
,
y
T
,
∀
j
=
1
,
…
,
4
,
which are calledsmoothed probabilities. We are going to show how these new probabilities can be computed from previous procedure (the same algorithm, although with some obvious changes, can be still used starting from procedure in Section 2.1).Since the last iteration of the algorithm gives us the probabilitiesP
(
S
T
=
j
∣
ψ
T
) for j
=
1
,
…
,
4, we can start from these values and use the following procedure by doing the two steps for every t
=
T
-
1
,
T
-
2
,
…
,
2,1.Step 1.
Fori
,
j
=
1
,
…
,
4, compute
(
*
)
P
S
t
=
i
,
S
t
+
1
=
j
∣
ψ
T
=
P
S
t
+
1
=
j
∣
ψ
T
P
S
t
=
i
∣
S
t
+
1
=
j
,
ψ
T
=
P
S
t
+
1
=
j
∣
ψ
T
P
S
t
=
i
∣
S
t
+
1
=
j
,
ψ
t
=
P
S
t
+
1
=
j
∣
ψ
T
P
S
t
=
i
,
S
t
+
1
=
j
∣
ψ
t
P
S
t
+
1
=
j
∣
ψ
t
.Remark 1.
Note that equality(
*
), that is,
(31)
P
S
t
=
i
∣
S
t
+
1
=
j
,
ψ
T
=
P
S
t
=
i
∣
S
t
+
1
=
j
,
ψ
t
,
holds only under a particular condition; namely,
(32)
f
h
t
+
1
,
T
∣
S
t
+
1
=
j
,
S
t
=
i
,
ψ
t
=
f
h
t
+
1
,
T
∣
S
t
+
1
=
j
,
ψ
t
,
where h
t
+
1
,
T
:
=
(
y
t
+
1
,
…
,
y
T
)
′ (see [6] for the proof). Equation (32) suggests that if S
t
+
1 were known, then y
t
+
1 would contain no information about S
t beyond that contained in S
t
+
1 and ψ
t and does not hold for every state-space model with regime switching (see, for example, [6, Ch.  5]) in which case the smoothing algorithm involves an approximation.Step 2.
Fori
=
1
,
…
,
4 compute
(33)
P
S
t
=
i
∣
ψ
T
=
∑
j
=
1
4
‍
P
S
t
=
i
,
S
t
+
1
=
j
∣
ψ
T
.


## 2.1. Serially Uncorrelated Data

Let us suppose thatS
=
S
t
t
∈
{
1
,
…
,
T
}, T
∈
N
+, is a discrete time stochastic process represented by a first-order Markov chain taking value in some set Ω
=
{
1
,
…
,
M
}, with transition probabilities:
(3)
p
i
j
∶
=
P
(
S
t
+
1
=
j
∣
S
t
=
i
,
S
t
-
1
=
i
t
-
1
,
…
,
S
0
=
i
0
)
=
P
(
S
t
+
1
=
j
∣
S
t
=
i
)
,
∀
t
∈
{
2
,
…
,
T
}
.
Then the state-space model we want to study is given by
(4)
y
t
=
μ
S
t
+
ϵ
t
,
t
=
1,2
,
…
,
T
,
ϵ
t
t
∈
{
1
,
…
,
T
}
are
i
.
i
.
d
.
N
(
0
,
σ
S
t
2
)
,
μ
S
t
=
μ
1
S
1
,
t
+
⋯
+
μ
M
S
M
,
t
,
σ
S
t
=
σ
1
S
1
,
t
+
⋯
+
σ
M
S
M
,
t
,
S
t
∈
{
1
,
…
,
M
}
with
probabilities
p
i
j
,
where the variables S
k
,
t, (
k
,
t
)
∈
{
1
,
…
,
M
}
×
{
1
,
…
,
T
} are introduced in order to have a slightly compact equation for μ
S
t and σ
S
t: in particular, S
k
,
t
=
1 if S
t
=
k, otherwise S
k
,
t
=
0, which implies that under regime k, for k
∈
{
1
,
…
,
M
}, parameters of mean and variance are given by μ
k and σ
k.Let us underline that in (4) the y
t are our observed data, for example, historical returns of a stock or some index time series, and we suppose that they arelocally distributed as Gaussian random variable in the sense that occasionally jumps could occur for both the mean μ
S
t and the variance σ
S
t
2. In particular, we assume that σ
1
<
σ
2
<
⋯
<
σ
M and we want to estimate these M unobserved values for standard deviation, as well as the M values for the switching mean. Note that we could also take μ as a constant, obtaining the so calledswitching variance problem or σ as a constant, having aswitching mean problem. The first one is in general more interesting in the analysis of financial time series, since variance is usually interpreted as an indicator for the market’s volatility.Given the model described by (4), the conditional density of y
t given S
t is Gaussian with mean μ
S
t and standard deviation σ
S
t; namely, its related probability density function reads as follow:
(5)
f
y
t
∣
S
t
(
x
)
=
1
2
π
σ
S
t
2
e
-
(
x
-
μ
S
t
)
2
/
2
σ
S
t
2
,
x
∈
R
,
and we are left with the problem of estimating both the expectations μ
1
,
…
,
μ
M and the standard deviations σ
1
,
…
,
σ
M parameters; a task that is standard to solve by maximizing the associated log-likelihood function ln
⁡
L defined as follows:
(6)
ln
⁡
L
:
=
∑
t
=
1
T
‍
ln
⁡
(
f
y
t
∣
S
t
)
.
A different and more realistic scenario is the one characterized by unobserved values for S
t. In such a case, it could be possible to consider the MSM-inference problem as a two-step procedure consisting in(i)
estimating the parameters of the model by maximizing the log-likelihood function,
(ii)
making inferences on the state variableS
t,  t
=
1
,
…
,
T conditional on the estimates obtained at previous point.Depending on the amount of information we can use inferencing onS
t, we have(i)
filtered probabilities that refer to inferences about S
t conditional on information up to time t, namely, with respect to ψ
t,
(ii)
smoothed probabilities that refer to inferences about S
t conditional to the whole sample (history), ψ
T.In what follows we describe a step-by-step algorithm which allows us to resolve the filtering problem for a sample of serially uncorrelated data. In particular we slightly generalize the approach given in [6], assuming that the state variable S
t belongs to a 4-state space set Ω
:
=
{
1,2
,
3,4
}, at every time t
=
1
,
…
,
T. Despite 2-state (expansion/contraction) and 3-state (low/medium/high volatility regime) models being the usual choices, we decided to consider 4-state MSM in order to refine the analysis with respect to volatility levels around the mean. A finer analysis can be also performed, even if one has to take into account the related nonlinear computational growth. We first define the log-likelihood function at time t as
(7)
l
(
θ
,
ψ
t
)
=
l
(
μ
1
,
…
,
μ
M
,
σ
1
,
…
,
σ
M
,
ψ
t
)
,
where θ
:
=
(
μ
1
,
…
,
μ
M
,
σ
1
,
…
,
σ
M
) is the vector of parameters that we want to estimate. Let us note that the θ is updated at every iteration, since we maximize with respect to the function l
(
θ
,
ψ
t
) at every stage of our step-by-step procedure. In particular the calibrating procedure reads as follows.Inputs.(i)
Putl
(
θ
)
:
=
l
(
θ
,
ψ
0
)
=
0.
(ii)
Compute thetransition probabilities of the homogeneous Markov chain underlying the state variable; that is,
(8)
p
i
j
=
P
(
S
t
=
j
∣
S
t
-
1
=
i
)
,
i
,
j
∈
Ω
.
Since in the applications we can only count on return time series, we first have to calibrate with respect to the transition probabilitiesp
i
j:
(1)
choose 4 valuesσ
1
2
<
σ
2
2
<
σ
3
2
<
σ
4
2 (e.g., σ
1
=
0.1, σ
2
=
0.2, σ
3
=
0.3, and σ
4
=
0.4) and a positive, arbitrary small, constant δ
>
0, for example, δ
=
0.01,
(2)
compute for everyj
=
1,2
,
3,4 and for every t
=
1
,
…
,
T the values
(9)
b
j
≔
Φ
y
t
+
δ
-
Φ
y
t
-
δ
=
∫
y
t
-
δ
y
t
+
δ
‍
1
2
π
σ
j
2
exp
⁡
x
-
μ
j
2
2
σ
j
2
,
(3)
simulate a value in{
1,2
,
3,4
} for S
t at each time from the discrete probability vector
(10)
b
1
∑
j
=
1
4
‍
b
j
,
b
2
∑
j
=
1
4
‍
b
j
,
b
3
∑
j
=
1
4
‍
b
j
,
b
4
∑
j
=
1
4
‍
b
j
,
(4)
set the transition probabilitiesp
i
j just by counting the number of transition from state i to state j, for i
,
j
=
1
,
…
,
4, in order to obtain the followingtransition matrix:
(11)
P
*
=
p
11
p
12
p
13
p
14
p
21
p
22
p
23
p
24
p
31
p
32
p
33
p
34
p
41
p
42
p
43
p
44
.
(iii)
Compute thesteady-state probabilities
(12)
π
(
0
)
=
P
S
0
=
1
∣
ψ
0
,
…
,
P
S
0
=
4
∣
ψ
0
.
Let us note that, by definition, ifπ
(
t
) is a 4
×
1 vector of steady-state probabilities, then π
(
t
+
1
)
=
π
(
t
) for every t
=
1
,
…
,
T; moreover π
(
t
+
1
)
=
P
*
π
(
t
), and (see for example, [6, pag.  71]) we also have that π
(
t
)
=
(
A
T
A
)
-
1
A
T
0
4
1, where 0
4 is a 4
×
1 matrix of zeros and A
=
I
d
4
-
P
*
1
4, I
d
4 is the four dimensional identity matrix, while 1
4
=
(
1,1
,
1,1
), that is, the vector of steady-state probabilities is the last column of the matrix (
A
T
A
)
-
1
A
T.Next we perform the following steps, fort
=
1
,
…
,
T.Step 1.
The probability ofS
t conditional to information set at time t
-
1 is given by
(13)
P
S
t
=
j
∣
ψ
t
-
1
=
∑
i
=
1
4
‍
p
i
j
P
S
t
-
1
=
i
∣
ψ
t
-
1
,
j
=
1
,
…
,
4
.Step 2.
Compute the joint density ofy
t and S
t conditional to the information set ψ
t
-
1
(14)
f
y
t
,
S
t
=
j
∣
ψ
t
-
1
=
f
y
t
∣
S
t
=
j
,
ψ
t
-
1
×
P
S
t
=
j
∣
ψ
t
-
1
,
j
=
1
,
…
,
4
.
The marginal density of y
t is given by the sum of the joint density over all values of S
t
(15)
f
y
t
∣
ψ
t
-
1
=
∑
i
=
1
4
‍
f
y
t
∣
S
t
=
i
,
ψ
t
-
1
P
S
t
=
i
∣
ψ
t
-
1
.Step 3.
Update the log-likelihood function at timet in the following way:
(16)
l
θ
,
ψ
t
=
l
θ
,
ψ
t
-
1
+
ln
⁡
f
y
t
∣
ψ
t
-
1
,
and maximize l
(
θ
,
ψ
t
) with respect to θ
=
(
μ
1
,
…
,
μ
4
,
σ
1
,
…
,
σ
4
), under the condition σ
1
<
σ
2
<
σ
3
<
σ
4, to find the maximum likelihood estimator θ
^ for the next time period.Step 4.
Oncey
t is observed at the end of the tth iteration, we canupdate the probability term:
(17)
P
S
t
=
j
∣
ψ
t
=
f
y
t
S
t
=
j
,
ψ
t
-
1
P
S
t
=
j
ψ
t
-
1
∑
i
=
1
4
‍
f
y
t
∣
S
t
=
i
,
ψ
t
-
1
P
S
t
=
i
∣
ψ
t
-
1
,
hhhhhhhhhhhhhhhhhhhhhhih
j
=
1
,
…
,
4
,
where both ψ
t
=
{
ψ
t
-
1
,
y
t
} and f
(
y
t
∣
S
t
=
i
,
ψ
t
-
1
) are computed with respect to the estimator θ
^
=
(
μ
^
1
,
…
,
μ
^
4
,
σ
^
1
,
…
,
σ
^
4
).


## 2.2. Serially Correlated Data

In some cases it is possible to argue and mathematically test by, for example, the Durbin-Watson statistics or Breusch-Godfrey test, for the presence of a serial correlation (orautocorrelation) between data belonging to a certain time series of interest. Such a characteristic is often analyzed in signal processing scenario, but examples can be also found in economic, meteorological, or sociological data sets, especially in connection with autocorrelation of errors in related forecasting procedure. In particular if we suppose that the observed variable y
t linearly depends on its previous value, then we obtain afirst-order autoregressive pattern and the following state-space model applies:
(18)
y
t
-
μ
S
t
=
ϕ
(
y
t
-
1
-
μ
S
t
-
1
)
+
ϵ
t
,
t
=
1,2
,
…
,
T
,
ϵ
t
t
∈
{
1
,
…
,
T
}
are
i
.
i
.
d
.
N
(
0
,
σ
S
t
2
)
,
μ
S
t
=
μ
1
S
1
,
t
+
⋯
+
μ
M
S
M
,
t
,
σ
S
t
=
σ
1
S
1
,
t
+
⋯
+
σ
M
S
M
,
t
,
S
t
∈
1
,
…
,
M
with
probabilities
p
i
j
,
where p
i
j
:
=
P
(
S
t
=
j
∣
S
t
-
1
=
i
), i
,
j
=
1
,
…
,
M and S
k
,
t, and (
k
,
t
)
∈
{
1
,
…
,
M
}
×
{
1
,
…
,
T
} are the same variables introduced in the previous section; that is, S
k
,
t
=
1 if S
t
=
k, otherwise S
k
,
t
=
0.In this situation, if the stateS
t is known for every t
=
1
,
…
,
T, we need S
t and S
t
-
1 to compute the density of y
t conditional to past information ψ
t
-
1, indeed we have
(19)
ln
⁡
L
:
=
∑
t
=
1
T
‍
ln
⁡
f
y
t
∣
ψ
t
-
1
,
S
t
,
S
t
-
1
,
where
(20)
f
y
t
∣
ψ
t
-
1
,
S
t
,
S
t
-
1
(
x
)
=
1
2
π
σ
S
t
2
e
-
x
-
μ
S
t
-
ϕ
y
t
-
1
-
μ
S
t
-
1
2
/
2
σ
S
t
2
,
x
∈
R
.IfS
t are unobserved (and as before we assume that the state variable can take the four values {
1,2
,
3,4
}), we apply the following algorithm in order to resolve the filtering problem for a sample of serially correlated data.Inputs.(i)
Putl
(
θ
)
:
=
l
(
θ
,
ψ
0
)
=
0.
(ii)
Compute the transition probabilities(21)
p
i
j
=
P
S
t
=
j
∣
S
t
-
1
=
i
,
i
,
j
=
1,2
,
3,4
.
We apply the same trick as before, but firstly we have to estimate the parameterϕ: in order to obtain this value we can use the least square methods (see for example, [7]) that is,
(22)
ϕ
^
:
=
∑
t
=
1
T
‍
y
t
y
t
+
1
∑
t
=
1
T
‍
y
t
2
.
Then we computez
t
:
=
y
t
-
ϕ
y
t
-
1 for every t
=
1
,
…
,
T and consider the values b
j
:
=
Φ
(
z
t
+
δ
)
-
Φ
(
z
t
-
δ
) (we apply the Normal distribution function to z
t
+
δ instead of y
t
+
δ, as done before).
(iii)
Compute the steady-state probabilities(23)
π
(
0
)
=
P
S
0
=
1
∣
ψ
0
,
…
,
P
S
0
=
4
∣
ψ
0
,
taking the last column of the matrix(
A
T
A
)
-
1
A
T (see procedure in Section 2.1 for details).Next perform the following steps for t
=
1
,
…
,
T.Step 1.
Compute the probabilities ofS
t conditional to information set at time t
-
1, for j
=
1
,
…
,
4
(24)
P
S
t
=
j
∣
ψ
t
-
1
=
∑
i
=
1
4
‍
P
S
t
=
j
,
S
t
-
1
=
i
∣
ψ
t
-
1
=
∑
i
=
1
4
‍
p
i
j
P
S
t
-
1
=
i
∣
ψ
t
-
1
.Step 2.
Compute the joint density ofy
t, S
t, and S
t
-
1 given ψ
t
-
1:
(25)
f
y
t
,
S
t
,
S
t
-
1
∣
ψ
t
-
1
=
f
y
t
∣
S
t
=
j
,
S
t
-
1
=
i
,
ψ
t
-
1
×
P
S
t
=
j
,
S
t
-
1
=
i
∣
ψ
t
-
1
,
where f
(
y
t
∣
S
t
=
j
,
S
t
-
1
=
i
,
ψ
t
-
1
) is given by (20) and P
(
S
t
=
j
,
S
t
-
1
=
i
∣
ψ
t
-
1
) is computed in Step 1. The marginal density of y
t conditional on ψ
t
-
1 is obtained by summing over all values of S
t and S
t
-
1:
(26)
f
y
t
∣
ψ
t
-
1
=
∑
j
=
1
4
‍
∑
i
=
1
4
‍
f
y
t
∣
S
t
=
j
,
S
t
-
1
=
i
,
ψ
t
-
1
hhhhl
×
P
S
t
=
j
,
S
t
-
1
=
i
∣
ψ
t
-
1
.Step 3.
The log-likelihood function at timet is again
(27)
l
θ
,
ψ
t
=
l
θ
,
ψ
t
-
1
+
ln
⁡
f
y
t
∣
ψ
t
-
1
,
and it can be maximized with respect to θ
=
(
μ
1
,
…
,
μ
4
,
σ
1
,
…
,
σ
4
), under condition σ
1
<
σ
2
<
σ
3
<
σ
4, giving the maximum likelihood estimator θ
^ for the next time period.Step 4.
Update the joint probabilities ofS
t and S
t
-
1 conditional to the new information set ψ
t, using the estimator θ
^ computed in Step 3 by maximizing the log-likelihood function l
(
θ
,
ψ
t
):
(28)
P
S
t
=
j
,
S
t
-
1
=
i
∣
ψ
t
=
f
y
t
,
S
t
=
j
,
S
t
-
1
=
i
∣
ψ
t
-
1
f
y
t
∣
ψ
t
-
1
,
i
,
j
=
1
,
…
,
4
.
Then compute the updated probabilities ofS
t given ψ
t by summing the joint probabilities over S
t
-
1 as follows:
(29)
P
S
t
=
j
∣
ψ
t
=
∑
i
=
1
4
‍
P
S
t
=
j
,
S
t
-
1
=
i
∣
ψ
t
,
∀
j
=
1
,
…
,
4
.The Smoothing Algorithm. Once we have run this procedure, we are provided with thefiltered probabilities, that is, the values P
(
S
t
=
j
∣
ψ
t
) for j
=
1
,
…
,
4 and for each t
=
1
,
…
,
T (in addition to the estimator θ
^).Sometimes it is required to estimate probabilities ofS
t given the whole sample information; that is,
(30)
P
S
t
=
j
∣
ψ
T
=
P
S
t
=
j
∣
y
1
,
…
,
y
T
,
∀
j
=
1
,
…
,
4
,
which are calledsmoothed probabilities. We are going to show how these new probabilities can be computed from previous procedure (the same algorithm, although with some obvious changes, can be still used starting from procedure in Section 2.1).Since the last iteration of the algorithm gives us the probabilitiesP
(
S
T
=
j
∣
ψ
T
) for j
=
1
,
…
,
4, we can start from these values and use the following procedure by doing the two steps for every t
=
T
-
1
,
T
-
2
,
…
,
2,1.Step 1.
Fori
,
j
=
1
,
…
,
4, compute
(
*
)
P
S
t
=
i
,
S
t
+
1
=
j
∣
ψ
T
=
P
S
t
+
1
=
j
∣
ψ
T
P
S
t
=
i
∣
S
t
+
1
=
j
,
ψ
T
=
P
S
t
+
1
=
j
∣
ψ
T
P
S
t
=
i
∣
S
t
+
1
=
j
,
ψ
t
=
P
S
t
+
1
=
j
∣
ψ
T
P
S
t
=
i
,
S
t
+
1
=
j
∣
ψ
t
P
S
t
+
1
=
j
∣
ψ
t
.Remark 1.
Note that equality(
*
), that is,
(31)
P
S
t
=
i
∣
S
t
+
1
=
j
,
ψ
T
=
P
S
t
=
i
∣
S
t
+
1
=
j
,
ψ
t
,
holds only under a particular condition; namely,
(32)
f
h
t
+
1
,
T
∣
S
t
+
1
=
j
,
S
t
=
i
,
ψ
t
=
f
h
t
+
1
,
T
∣
S
t
+
1
=
j
,
ψ
t
,
where h
t
+
1
,
T
:
=
(
y
t
+
1
,
…
,
y
T
)
′ (see [6] for the proof). Equation (32) suggests that if S
t
+
1 were known, then y
t
+
1 would contain no information about S
t beyond that contained in S
t
+
1 and ψ
t and does not hold for every state-space model with regime switching (see, for example, [6, Ch.  5]) in which case the smoothing algorithm involves an approximation.Step 2.
Fori
=
1
,
…
,
4 compute
(33)
P
S
t
=
i
∣
ψ
T
=
∑
j
=
1
4
‍
P
S
t
=
i
,
S
t
+
1
=
j
∣
ψ
T
.


## 3. The Gibbs Sampling Approach

### 3.1. An Introduction to Bayesian Inference

Under the general titleBayesian inference we can collect a large number of different concrete procedures; nevertheless they are all based on smart use of the Bayes’ rule which is used to update the probability estimate for a hypothesis as additional evidence is learned (see, for example, [8, 9]). In particular, within the Bayesian framework, the parameters, for example, let us collect them in a vector called θ, which characterize a certain statistic model and are treated as random variables with their own probability distributions; let us say f
(
θ
), which plays the role of aprior distribution since it is defined before taking into account the sample data y
~. Therefore, exploiting the Bayes’ theorem and denoting by f
(
y
~
∣
θ
) the likelihood of y
~ of the interested statistic model, we have that
(34)
f
θ
∣
y
~
=
f
y
~
∣
θ
f
θ
f
y
~
,
where f
(
θ
∣
y
~
) is thejoint posterior distribution of the parameters. The denominator f
(
y
~
) defines the marginal likelihood of y
~ and can be taken as a constant, obtaining the proportion
(35)
f
θ
∣
y
~
∝
f
y
~
∣
θ
f
(
θ
)
.
It is straightforward to note that the most critical part of the Bayesian inference procedure relies in the choice of a suitable prior distribution, since it has to agree with parameters constraints. An effective answer to latter issue is given by the so calledconjugate prior distribution, namely the distribution obtained when the conjugate prior is combined with the likelihood function. Let us note that the posterior distribution f
(
θ
∣
y
~
) is in the same family as the prior distribution.As an example, if the likelihood function is Gaussian, it can be shown that the conjugate prior for the meanμ is the Gaussian distribution, whereas the conjugate prior for the variance is the inverted Gamma distribution (see, for example, [9, 10]).


### 3.2. Gibbs Sampling

A general problem in Statistics concerns the question of how a sequence of observations which cannot be directly sampled, can be simulated, for example, by mean of some multivariate probability distribution, with a prefixed precision degree of accuracy. Such kind of problems can be successfully attacked by Monte Carlo Markov Chain (MCMC) simulation methods, see, for example, [11–13], and in particular using the so calledGibbs Sampling technique which allows to approximate joint and marginal distributions by sampling from conditional distributions, see, for example, [14–16].Let us suppose that we have the joint density ofk random variables, for example, f
=
f
(
z
1
,
z
2
,
…
,
z
k
), fix t
∈
{
1
,
…
,
k
} and that we are interested in in obtaining characteristics of the z
t-marginal, namely
(36)
f
z
t
=
∫
‍
⋯
∫
‍
f
(
z
1
,
z
2
,
…
,
z
k
)
d
z
1
⋯
d
z
t
-
1
d
z
t
+
1
⋯
d
z
k
,
such as the related mean and/or variance. In those cases when the joint density is not given, or the above integral turns out to be difficult to treat, for example, an explicit solution does not exist, but we know the complete set of conditional densities, denoted by f
(
z
t
∣
z
j
≠
t
), t
=
1,2
,
…
,
k, with z
j
≠
t
=
{
z
1
,
…
,
z
t
-
1
,
z
t
+
1
,
…
,
z
k
}, then the Gibbs Sampling method allows us to generate a sample z
1
j
,
z
2
j
,
…
,
z
k
j from the joint density f
(
z
1
,
z
2
,
…
,
z
k
) without requiring that we know either the joint density or the marginal densities. With the following procedure we recall the basic ideas on which the Gibbs Sampling approach is based given an arbitrary starting set of values (
z
2
0
,
…
,
z
k
0
).Step 1.
Drawz
1
1 from f
(
z
1
∣
z
2
0
,
…
,
z
k
0
).Step 2.
Drawz
2
1 from f
(
z
2
∣
z
1
1
,
z
3
0
,
…
,
z
k
0
).Step 3.
Drawz
3
1 from f
(
z
3
∣
z
1
1
,
z
2
1
,
z
4
0
,
…
,
z
k
0
).⋮Step k. Finally draw z
k
1 from f
(
z
k
∣
z
1
1
,
…
,
z
k
-
1
1
) to complete the first iteration.The steps from1 through k can be iterated J times to get (
z
1
j
,
z
2
j
,
…
,
z
k
j
), j
=
1,2
,
…
,
J.In [17] S. Geman and D. Geman showed that both the joint and marginal distributions of generated (
z
1
j
,
z
2
j
,
…
,
z
k
j
), converge, at an exponential rate, to the joint and marginal distributions of z
1
,
z
2
,
…
,
z
k, as J
→
∞. Thus the joint and marginal distributions of z
1
,
z
2
,
…
,
z
k can be approximated by the empirical distributions of M simulated values (
z
1
j
,
z
2
j
,
…
,
z
k
j
), j
=
L
+
1
,
…
,
L
+
M, where L is large enough to assure the convergence of the Gibbs sampler. Moreover M can be chosen to reach the required precision with respect to the empirical distribution of interest.In the MSM framework we do not have conditional distributionsf
(
z
t
∣
z
j
≠
t
),  t
=
1,2
,
…
,
k, and we are left with the problem of estimate parameters z
i, i
=
1
,
…
,
k. Latter problem can be solved exploiting Bayesian inference results, as we shall state in the next section.


### 3.3. Gibbs Sampling for Markov Switching Models

A major problem when dealing with inferences from Markov switching models relies in the fact that some parameters of the model are dependent on an unobserved variable, let us sayS
t. We saw that in the classical framework, inference on Markov switching models consists first in estimating the model’s unknown parameters via maximum likelihood, then inference on the unobserved Markov switching variable S
~
T
=
(
S
1
,
S
2
,
…
,
S
T
), conditional on the parameter estimates, has to be perfomed.In the Bayesian analysis, both the parameters of the model and the switching variablesS
t,  t
=
1
,
…
,
T are treated as random variables. Thus, inference on S
~
T is based on a joint distribution, no more on a conditional one. By employing Gibbs sampling techniques, Albert and Chib (see [14]) provided an easy to implement algorithm for the Bayesian analysis of Markov switching models. In particular in their work the parameters of the model and S
t,  t
=
1
,
…
,
T are treated as missing data, and they are generated from appropriate conditional distributions using Gibbs sampling method. As an example, let us consider the following simple model with two-state Markov switching mean and variance:
(37)
y
t
=
μ
S
t
+
ϵ
t
,
t
=
1,2
,
…
,
T
,
ϵ
t
t
∈
{
1
,
…
,
T
}
are
i
.
i
.
d
.
N
(
0
,
σ
S
t
2
)
,
μ
S
t
=
μ
0
+
μ
1
S
t
,
σ
S
t
=
σ
0
2
(
1
-
S
t
)
+
σ
1
2
S
t
=
σ
0
2
(
1
+
h
1
S
t
)
,
h
1
>
0
,
where S
t
∈
{
0,1
} with transition probabilities p
:
=
P
(
S
t
=
0
∣
S
t
-
1
=
0
), q
:
=
P
(
S
t
=
1
∣
S
t
-
1
=
1
). The Bayesian method consider both S
t, t
=
1
,
…
,
T and the model’s unknown parameters μ
0, μ
1, σ
0, σ
1, p and q, as random variables. In order to make inference about these T
+
6 variables, we need to derive the joint posterior density f
(
S
~
T
,
μ
0
,
μ
1
,
σ
0
2
,
σ
1
2
,
p
,
q
∣
ψ
T
), where ψ
T
:
=
(
y
1
,
y
2
,
…
,
y
T
) and S
~
T
=
(
S
1
,
S
2
,
…
,
S
T
). Namely the realization of the Gibbs sampling relies on the derivation of the distributions of each of the above T
+
6 variables, conditional on all the other variables. Therefore we can approximate the joint posterior density written above by running the following procedure L
+
M times, where L is an integer large enough to guarantee the desired convergence. Hence we have the following scheme.Step 1.
We can derive the distribution ofS
t, t
=
1
,
…
,
T conditional on the other parameters in two different ways.(1)
Single-move gibbs sampling: generate each S
t from f
(
S
t
∣
S
~
≠
t
,
μ
0
,
μ
1
,
σ
0
2
,
σ
1
2
,
p
,
q
,
ψ
T
), t
=
1
,
…
,
T, where S
~
≠
t
=
(
S
1
,
…
,
S
t
-
1
,
S
t
+
1
,
…
,
S
T
).
(2)
Multi-move gibbs sampling: generate the whole block S
~
T from f
(
S
~
T
∣
μ
0
,
μ
1
,
σ
0
2
,
σ
1
2
,
p
,
q
,
ψ
T
).Step 2.
Generate the transition probabilitiesp and q from f
(
p
,
q
∣
S
~
T
). Note that this distribution is conditioned only on S
~
T because we assume that p and q are independent of both the other parameters of the model and the data, ψ
T.
If we choose theBeta distribution as prior distribution for both p and q, we have that posterior distribution f
(
p
,
q
∣
S
~
T
)
=
f
(
p
,
q
)
L
(
p
,
q
∣
S
~
T
) is again a Beta distribution. So, Beta distribution is a conjugate prior for the likelihood of transition probabilities.Step 3.
Generateμ
0 and μ
1 from f
(
μ
0
,
μ
1
∣
S
~
T
,
σ
0
2
,
σ
1
2
,
p
,
q
,
ψ
T
). In this case the conjugate prior is theNormal distribution.Step 4.
Generateσ
0
2 and σ
1
2 from f
(
σ
0
2
,
σ
1
2
∣
S
~
T
,
μ
0
,
μ
1
,
p
,
q
,
ψ
T
). From definition of the model we have that σ
1
2
=
σ
0
2
(
1
+
h
1
): we can first generate σ
0
2 conditional on h
1, and then generate h
-
1
=
1
+
h
1 conditional on σ
0
2. We use in both cases theInverted Gamma distribution as conjugate prior for the parameters.For a more detailed description of these steps (see [6, pp.  211–218]). Here we examine only the so calledMulti-move Gibbs sampling, originally motivated by Carter and Kohn (see [15]) in the context of state space models and then implemented in [6] for a MSM. For the sake of simplicity, let us suppress the conditioning on model’s parameters and denote
(38)
f
S
~
T
∣
ψ
T
:
=
f
S
~
T
∣
μ
0
,
μ
1
,
σ
0
2
,
σ
1
2
,
p
,
q
,
ψ
T
.Using the Markov property of{
S
t
}
t
∈
{
1
,
…
,
T
} it can be seen that
(39)
f
S
~
T
∣
ψ
T
=
f
S
T
∣
ψ
T
∏
t
=
1
T
-
1
‍
f
S
t
∣
S
t
+
1
,
ψ
t
,
where f
(
S
T
∣
ψ
T
)
=
P
(
S
T
∣
ψ
T
) is provided by the last iteration of filtering algorithm (see Sections 2.1 and 2.2). Note that (39) suggests that we can first generate S
T conditional on ψ
T and then, for t
=
T
-
1
,
T
-
2
,
…
,
1, we can generate S
t conditional on ψ
t and S
t
+
1, namely we can run the following steps.Step 1.
Run the basic filter procedure to getf
(
S
t
∣
ψ
t
), t
=
1,2
,
…
,
T and save them; the last iteration of the filter gives us the probability distribution f
(
S
T
∣
ψ
T
), from which S
T is generated.Step 2.
Note that(40)
f
(
S
t
∣
S
t
+
1
,
ψ
t
)
=
f
(
S
t
,
S
t
+
1
∣
ψ
t
)
f
(
S
t
+
1
∣
ψ
t
)
=
f
(
S
t
+
1
∣
S
t
,
ψ
t
)
f
(
S
t
∣
ψ
t
)
f
(
S
t
+
1
∣
ψ
t
)
=
f
(
S
t
+
1
∣
S
t
)
f
(
S
t
∣
ψ
t
)
f
(
S
t
+
1
∣
ψ
t
)
∝
f
(
S
t
+
1
∣
S
t
)
f
(
S
t
∣
ψ
t
)
,
where f
(
S
t
+
1
∣
S
t
) is the transition probability and f
(
S
t
∣
ψ
t
) has been saved from Step 1. So we can generate S
t in the following way, first calculate
(41)
P
(
S
t
=
1
∣
S
t
+
1
,
ψ
t
)
=
f
(
S
t
+
1
∣
S
t
=
1
)
f
(
S
t
=
1
∣
ψ
t
)
∑
j
=
0
1
‍
f
(
S
t
+
1
∣
S
t
=
j
)
f
(
S
t
=
j
∣
ψ
t
)
,
and then generate S
t using a uniform distribution. For example, we generate a random number from a uniform distribution between 0 and 1; if this number is less than or equal to the calculated value of P
(
S
t
=
1
∣
S
t
+
1
,
ψ
t
), we set S
t
=
1, otherwise, S
t is set equal to 0.In view of applications, let us now consider the followingfour state MSM:
(42)
y
t
~
N
0
,
σ
S
t
2
,
t
=
1,2
,
…
,
T
,
σ
S
t
=
σ
1
S
1
,
t
+
σ
2
S
2
,
t
+
σ
3
S
3
,
t
+
σ
4
S
4
,
t
,
l
S
t
∈
{
1,2
,
3,4
}
with
transition
probabilities
p
i
j
:
=
P
(
S
t
=
j
∣
S
t
-
1
=
i
)
,
i
,
j
=
1,2
,
3,4
,
where S
k
,
t
=
1 if S
t
=
k, otherwise S
k
,
t
=
0. Note that this is a particular case of the model analysed in Section 2.1, where μ
t
=
0  ∀
t, hence we can perform the procedure referred to serially uncorrelated data taking μ
S
t
=
μ
=
0 to start the Gibbs sampling algorithm, therefore we haveStep 1.
GenerateS
~
T
:
=
(
S
1
,
S
2
,
…
,
S
T
) conditional on σ
~
2
:
=
(
σ
1
2
,
σ
2
2
,
σ
3
2
,
σ
4
2
),
(43)
p
~
:
=
p
11
,
p
12
,
p
13
,
p
21
,
p
22
,
p
23
,
p
31
,
p
32
,
p
33
,
p
41
,
p
42
,
p
43
,
ψ
T
:
=
y
1
,
y
2
,
…
,
y
T
.
For this purpose, we employ the Multi-move Gibbs sampling algorithm:(1)
run procedure in Section2.1 with μ
S
t
=
0 in order to get, from last iteration, f
(
S
T
∣
ψ
T
)
=
P
(
S
T
∣
ψ
T
),
(2)
recalling thatf
(
S
t
∣
S
t
+
1
,
ψ
t
)
∝
f
(
S
t
+
1
∣
S
t
)
f
(
S
t
∣
ψ
t
) for t
=
T
-
1
,
…
,
1, we can generate S
t from the vector of probabilities
(44)
P
(
S
t
=
1
∣
S
t
+
1
,
ψ
t
)
,
P
(
S
t
=
2
∣
S
t
+
1
,
ψ
t
)
,
h
P
(
S
t
=
3
∣
S
t
+
1
,
ψ
t
)
,
P
(
S
t
=
4
∣
S
t
+
1
,
ψ
t
)
,where, for i
=
1
,
…
,
4,
(45)
P
S
t
=
i
S
t
+
1
,
ψ
t
=
f
S
t
+
1
S
t
=
1
f
S
t
=
i
ψ
t
∑
j
=
1
3
‍
f
S
t
+
1
S
t
=
j
f
S
t
=
j
ψ
t
.Step 2.
Generateσ
~
2, conditional on S
~
T and the data ψ
T.
We want to impose the constraintσ
1
2
<
σ
2
2
<
σ
3
2
<
σ
4
2, so we redefine σ
S
t
2 in this way:
(46)
σ
S
t
2
=
σ
1
2
1
+
S
2
,
t
h
2
1
+
S
3
,
t
h
2
1
+
S
3
,
t
h
3
×
(
1
+
S
4
,
t
h
2
)
(
1
+
S
4
,
t
h
3
)
(
1
+
S
4
,
t
h
4
)
,
where h
j
>
0 for j
=
1
,
…
,
4, so that σ
2
2
:
=
σ
1
2
(
1
+
h
2
), σ
3
2
:
=
σ
1
2
(
1
+
h
2
)
(
1
+
h
3
) and σ
4
2
:
=
σ
1
2
(
1
+
h
2
)
(
1
+
h
3
)
(
1
+
h
4
). With this specification, we first generate σ
1
2, then generate h
-
2
=
1
+
h
2, h
-
3
=
1
+
h
3 and h
-
4
=
1
+
h
4 to obtain σ
2
2, σ
3
2 and σ
4
2 indirectly.Generating
σ
1
2
, Conditional on
h
2
,
h
3
and
h
4. Define for t
=
1
,
…
,
T
(47)
Y
t
1
:
=
y
t
(
1
+
S
2
,
t
h
2
)
(
1
+
S
3
,
t
h
2
)
(
1
+
S
3
,
t
h
3
)
(
1
+
S
4
,
t
h
2
)
(
1
+
S
4
,
t
h
3
)
(
1
+
S
4
,
t
h
4
)and takeY
1
t
~
N
(
0
,
σ
1
2
), in (42). By choosing aninverted Gamma prior distribution, that is, f
(
σ
1
2
∣
h
2
,
h
3
,
h
4
)
~
IG
(
ν
1
/
2
,
δ
1
/
2
) where ν
1
,
δ
1 are the known prior hyperparameters, it can be shown that the conditional posterior distribution from which we generate σ
1
2 is given by:
(48)
σ
1
2
∣
ψ
T
,
S
~
T
,
h
2
,
h
3
,
h
4
~
IG
ν
1
+
T
2
,
δ
1
+
∑
t
=
1
T
‍
Y
t
1
2
.Generating
h
-
2
Conditional on
σ
1
2
,
h
3
and
h
4. Note that the likelihood function of h
2 depends only on the values of y
t for which S
t
∈
{
2,3
,
4
}. Therefore, take y
t
(
1
)
:
=
{
y
t
∣
S
t
∈
{
2,3
,
4
}
,
t
=
1
,
…
,
T
} and denote with T
2 the size of this sample. Then define
(49)
Y
t
2
:
=
y
t
(
1
)
σ
1
2
(
1
+
S
3
,
t
h
3
)
(
1
+
S
4
,
t
h
3
)
(
1
+
S
4
,
t
h
4
)
,
hence, for the observation in which S
t
=
2, 3 or 4, we have Y
2
t
~
N
(
0
,
h
-
2
). If we choose an inverted Gamma distribution with parameters ν
2
,
δ
2 for the prior, we obtain h
-
2
=
1
+
h
2 from the following posterior distribution:
(50)
h
-
2
∣
ψ
T
,
S
~
T
,
σ
1
2
,
h
3
,
h
4
~
IG
ν
2
+
T
2
2
,
δ
2
+
∑
t
=
1
T
2
‍
Y
t
2
2
.
In case h
-
2
>
1 put h
2
=
h
-
2
-
1 and σ
2
2
:
=
σ
1
2
(
1
+
h
2
). Otherwise reiterate this step.Generating
h
-
3
Conditional on
σ
1
2
,
h
2
and
h
4. Operate in a similar way as above. In particular if we define y
t
(
2
)
:
=
{
y
t
∣
S
t
∈
{
3,4
}
,
t
=
1
,
…
,
T
} we will obtain
(51)
Y
t
3
:
=
y
t
(
2
)
σ
1
2
(
1
+
S
3
,
t
h
2
)
(
1
+
S
4
,
t
h
2
)
(
1
+
S
4
,
t
h
4
)
~
N
(
0
,
h
-
3
)
.Generating
h
-
4
Conditional on
σ
1
2
,
h
2
and
h
3. Operate in a similar way as above. In particular if we define y
t
(
3
)
:
=
{
y
t
∣
S
t
=
4
,
t
=
1
,
…
,
T
} we will have
(52)
Y
t
4
:
=
y
t
(
3
)
σ
1
2
(
1
+
S
4
,
t
h
2
)
(
1
+
S
4
,
t
h
3
)
~
N
(
0
,
h
-
4
)
.Step 3.
Generatep
~ conditional on S
~
T. In order to generate the transition probabilities we exploit the properties of theprior Beta distribution. Let us first define
(53)
p
-
i
i
:
=
P
(
S
t
≠
i
∣
S
t
-
1
=
i
)
=
1
-
p
i
i
,
i
=
1,2
,
3,4
;
p
-
i
j
:
=
P
(
S
t
=
j
∣
S
t
-
1
=
i
,
S
t
≠
i
)
,
i
≠
j
.
Hence we have that
(54)
p
i
j
:
=
P
(
S
t
=
j
∣
S
t
-
1
=
i
)
=
P
(
S
t
=
j
∣
S
t
-
1
=
i
,
S
t
≠
i
)
P
(
S
t
≠
i
∣
S
t
-
1
=
i
)
=
p
-
i
j
(
1
-
p
i
i
)
∀
i
≠
j
.
Given S
~
T, let n
i
j, i
,
j
=
1,2
,
3,4 be the total number of transitions from state S
t
-
1
=
i to S
t
=
j,  t
=
2,3
,
…
,
T and n
-
i
j the number of transitions from state S
t
-
1
=
i to S
t
≠
j.
Begin with the generation of probabilitiesp
i
i,  i
=
1,2
,
3,4 by taking the Beta distribution as conjugate prior: if we take p
i
i
~
Beta
(
u
i
i
,
u
-
i
i
), where u
i
i and u
-
i
i are the known hyperparameters of the priors, the posterior distribution of p
i
i given S
~
T still belongs to the Beta family distributions, that is,
(55)
p
i
i
∣
S
~
T
~
Beta
(
u
i
i
+
n
i
i
,
u
-
i
i
+
n
-
i
i
)
,
i
=
1,2
,
3,4
.
The others parameters, that is, p
i
j for i
≠
j and j
=
1,2
,
3, can be computed from the above equation p
i
j
=
p
-
i
j
(
1
-
p
i
i
), where p
-
i
j are generated from the following posterior Beta distribution:
(56)
p
-
i
j
∣
S
~
T
~
Beta
(
u
i
j
+
n
i
j
,
u
-
i
j
+
n
-
i
j
)
.
For example, given that p
11 is generated, we can obtain p
12 and p
13 by considering
(57)
p
-
12
∣
S
~
T
~
Beta
(
u
12
+
n
12
,
u
-
12
+
n
-
12
)
,
p
-
13
∣
S
~
T
~
Beta
(
u
13
+
n
13
,
u
-
13
+
n
-
13
)
,
where n
-
12
=
n
13
+
n
14 and n
-
13
=
n
12
+
n
14. Finally, given p
11, p
12 and p
13 generated in this way, we have p
14
=
1
-
p
11
-
p
12
-
p
13.Remark 2.
When we do not have any information about priors distribution we employ hyperparametersu
i
j
=
0.5, i
,
j
=
1,2
,
3,4. Usually we know that elements of the matrix diagonal in the transition matrix are bigger than elements out of the diagonal, because in a financial framework regime switching happens only occasionally: in this case, since we want p
i
i close to 1 and p
i
j, i
≠
j, close to 0, we will choose u
i
i bigger than u
-
i
i.


## 3.1. An Introduction to Bayesian Inference

Under the general titleBayesian inference we can collect a large number of different concrete procedures; nevertheless they are all based on smart use of the Bayes’ rule which is used to update the probability estimate for a hypothesis as additional evidence is learned (see, for example, [8, 9]). In particular, within the Bayesian framework, the parameters, for example, let us collect them in a vector called θ, which characterize a certain statistic model and are treated as random variables with their own probability distributions; let us say f
(
θ
), which plays the role of aprior distribution since it is defined before taking into account the sample data y
~. Therefore, exploiting the Bayes’ theorem and denoting by f
(
y
~
∣
θ
) the likelihood of y
~ of the interested statistic model, we have that
(34)
f
θ
∣
y
~
=
f
y
~
∣
θ
f
θ
f
y
~
,
where f
(
θ
∣
y
~
) is thejoint posterior distribution of the parameters. The denominator f
(
y
~
) defines the marginal likelihood of y
~ and can be taken as a constant, obtaining the proportion
(35)
f
θ
∣
y
~
∝
f
y
~
∣
θ
f
(
θ
)
.
It is straightforward to note that the most critical part of the Bayesian inference procedure relies in the choice of a suitable prior distribution, since it has to agree with parameters constraints. An effective answer to latter issue is given by the so calledconjugate prior distribution, namely the distribution obtained when the conjugate prior is combined with the likelihood function. Let us note that the posterior distribution f
(
θ
∣
y
~
) is in the same family as the prior distribution.As an example, if the likelihood function is Gaussian, it can be shown that the conjugate prior for the meanμ is the Gaussian distribution, whereas the conjugate prior for the variance is the inverted Gamma distribution (see, for example, [9, 10]).


## 3.2. Gibbs Sampling

A general problem in Statistics concerns the question of how a sequence of observations which cannot be directly sampled, can be simulated, for example, by mean of some multivariate probability distribution, with a prefixed precision degree of accuracy. Such kind of problems can be successfully attacked by Monte Carlo Markov Chain (MCMC) simulation methods, see, for example, [11–13], and in particular using the so calledGibbs Sampling technique which allows to approximate joint and marginal distributions by sampling from conditional distributions, see, for example, [14–16].Let us suppose that we have the joint density ofk random variables, for example, f
=
f
(
z
1
,
z
2
,
…
,
z
k
), fix t
∈
{
1
,
…
,
k
} and that we are interested in in obtaining characteristics of the z
t-marginal, namely
(36)
f
z
t
=
∫
‍
⋯
∫
‍
f
(
z
1
,
z
2
,
…
,
z
k
)
d
z
1
⋯
d
z
t
-
1
d
z
t
+
1
⋯
d
z
k
,
such as the related mean and/or variance. In those cases when the joint density is not given, or the above integral turns out to be difficult to treat, for example, an explicit solution does not exist, but we know the complete set of conditional densities, denoted by f
(
z
t
∣
z
j
≠
t
), t
=
1,2
,
…
,
k, with z
j
≠
t
=
{
z
1
,
…
,
z
t
-
1
,
z
t
+
1
,
…
,
z
k
}, then the Gibbs Sampling method allows us to generate a sample z
1
j
,
z
2
j
,
…
,
z
k
j from the joint density f
(
z
1
,
z
2
,
…
,
z
k
) without requiring that we know either the joint density or the marginal densities. With the following procedure we recall the basic ideas on which the Gibbs Sampling approach is based given an arbitrary starting set of values (
z
2
0
,
…
,
z
k
0
).Step 1.
Drawz
1
1 from f
(
z
1
∣
z
2
0
,
…
,
z
k
0
).Step 2.
Drawz
2
1 from f
(
z
2
∣
z
1
1
,
z
3
0
,
…
,
z
k
0
).Step 3.
Drawz
3
1 from f
(
z
3
∣
z
1
1
,
z
2
1
,
z
4
0
,
…
,
z
k
0
).⋮Step k. Finally draw z
k
1 from f
(
z
k
∣
z
1
1
,
…
,
z
k
-
1
1
) to complete the first iteration.The steps from1 through k can be iterated J times to get (
z
1
j
,
z
2
j
,
…
,
z
k
j
), j
=
1,2
,
…
,
J.In [17] S. Geman and D. Geman showed that both the joint and marginal distributions of generated (
z
1
j
,
z
2
j
,
…
,
z
k
j
), converge, at an exponential rate, to the joint and marginal distributions of z
1
,
z
2
,
…
,
z
k, as J
→
∞. Thus the joint and marginal distributions of z
1
,
z
2
,
…
,
z
k can be approximated by the empirical distributions of M simulated values (
z
1
j
,
z
2
j
,
…
,
z
k
j
), j
=
L
+
1
,
…
,
L
+
M, where L is large enough to assure the convergence of the Gibbs sampler. Moreover M can be chosen to reach the required precision with respect to the empirical distribution of interest.In the MSM framework we do not have conditional distributionsf
(
z
t
∣
z
j
≠
t
),  t
=
1,2
,
…
,
k, and we are left with the problem of estimate parameters z
i, i
=
1
,
…
,
k. Latter problem can be solved exploiting Bayesian inference results, as we shall state in the next section.


## 3.3. Gibbs Sampling for Markov Switching Models

A major problem when dealing with inferences from Markov switching models relies in the fact that some parameters of the model are dependent on an unobserved variable, let us sayS
t. We saw that in the classical framework, inference on Markov switching models consists first in estimating the model’s unknown parameters via maximum likelihood, then inference on the unobserved Markov switching variable S
~
T
=
(
S
1
,
S
2
,
…
,
S
T
), conditional on the parameter estimates, has to be perfomed.In the Bayesian analysis, both the parameters of the model and the switching variablesS
t,  t
=
1
,
…
,
T are treated as random variables. Thus, inference on S
~
T is based on a joint distribution, no more on a conditional one. By employing Gibbs sampling techniques, Albert and Chib (see [14]) provided an easy to implement algorithm for the Bayesian analysis of Markov switching models. In particular in their work the parameters of the model and S
t,  t
=
1
,
…
,
T are treated as missing data, and they are generated from appropriate conditional distributions using Gibbs sampling method. As an example, let us consider the following simple model with two-state Markov switching mean and variance:
(37)
y
t
=
μ
S
t
+
ϵ
t
,
t
=
1,2
,
…
,
T
,
ϵ
t
t
∈
{
1
,
…
,
T
}
are
i
.
i
.
d
.
N
(
0
,
σ
S
t
2
)
,
μ
S
t
=
μ
0
+
μ
1
S
t
,
σ
S
t
=
σ
0
2
(
1
-
S
t
)
+
σ
1
2
S
t
=
σ
0
2
(
1
+
h
1
S
t
)
,
h
1
>
0
,
where S
t
∈
{
0,1
} with transition probabilities p
:
=
P
(
S
t
=
0
∣
S
t
-
1
=
0
), q
:
=
P
(
S
t
=
1
∣
S
t
-
1
=
1
). The Bayesian method consider both S
t, t
=
1
,
…
,
T and the model’s unknown parameters μ
0, μ
1, σ
0, σ
1, p and q, as random variables. In order to make inference about these T
+
6 variables, we need to derive the joint posterior density f
(
S
~
T
,
μ
0
,
μ
1
,
σ
0
2
,
σ
1
2
,
p
,
q
∣
ψ
T
), where ψ
T
:
=
(
y
1
,
y
2
,
…
,
y
T
) and S
~
T
=
(
S
1
,
S
2
,
…
,
S
T
). Namely the realization of the Gibbs sampling relies on the derivation of the distributions of each of the above T
+
6 variables, conditional on all the other variables. Therefore we can approximate the joint posterior density written above by running the following procedure L
+
M times, where L is an integer large enough to guarantee the desired convergence. Hence we have the following scheme.Step 1.
We can derive the distribution ofS
t, t
=
1
,
…
,
T conditional on the other parameters in two different ways.(1)
Single-move gibbs sampling: generate each S
t from f
(
S
t
∣
S
~
≠
t
,
μ
0
,
μ
1
,
σ
0
2
,
σ
1
2
,
p
,
q
,
ψ
T
), t
=
1
,
…
,
T, where S
~
≠
t
=
(
S
1
,
…
,
S
t
-
1
,
S
t
+
1
,
…
,
S
T
).
(2)
Multi-move gibbs sampling: generate the whole block S
~
T from f
(
S
~
T
∣
μ
0
,
μ
1
,
σ
0
2
,
σ
1
2
,
p
,
q
,
ψ
T
).Step 2.
Generate the transition probabilitiesp and q from f
(
p
,
q
∣
S
~
T
). Note that this distribution is conditioned only on S
~
T because we assume that p and q are independent of both the other parameters of the model and the data, ψ
T.
If we choose theBeta distribution as prior distribution for both p and q, we have that posterior distribution f
(
p
,
q
∣
S
~
T
)
=
f
(
p
,
q
)
L
(
p
,
q
∣
S
~
T
) is again a Beta distribution. So, Beta distribution is a conjugate prior for the likelihood of transition probabilities.Step 3.
Generateμ
0 and μ
1 from f
(
μ
0
,
μ
1
∣
S
~
T
,
σ
0
2
,
σ
1
2
,
p
,
q
,
ψ
T
). In this case the conjugate prior is theNormal distribution.Step 4.
Generateσ
0
2 and σ
1
2 from f
(
σ
0
2
,
σ
1
2
∣
S
~
T
,
μ
0
,
μ
1
,
p
,
q
,
ψ
T
). From definition of the model we have that σ
1
2
=
σ
0
2
(
1
+
h
1
): we can first generate σ
0
2 conditional on h
1, and then generate h
-
1
=
1
+
h
1 conditional on σ
0
2. We use in both cases theInverted Gamma distribution as conjugate prior for the parameters.For a more detailed description of these steps (see [6, pp.  211–218]). Here we examine only the so calledMulti-move Gibbs sampling, originally motivated by Carter and Kohn (see [15]) in the context of state space models and then implemented in [6] for a MSM. For the sake of simplicity, let us suppress the conditioning on model’s parameters and denote
(38)
f
S
~
T
∣
ψ
T
:
=
f
S
~
T
∣
μ
0
,
μ
1
,
σ
0
2
,
σ
1
2
,
p
,
q
,
ψ
T
.Using the Markov property of{
S
t
}
t
∈
{
1
,
…
,
T
} it can be seen that
(39)
f
S
~
T
∣
ψ
T
=
f
S
T
∣
ψ
T
∏
t
=
1
T
-
1
‍
f
S
t
∣
S
t
+
1
,
ψ
t
,
where f
(
S
T
∣
ψ
T
)
=
P
(
S
T
∣
ψ
T
) is provided by the last iteration of filtering algorithm (see Sections 2.1 and 2.2). Note that (39) suggests that we can first generate S
T conditional on ψ
T and then, for t
=
T
-
1
,
T
-
2
,
…
,
1, we can generate S
t conditional on ψ
t and S
t
+
1, namely we can run the following steps.Step 1.
Run the basic filter procedure to getf
(
S
t
∣
ψ
t
), t
=
1,2
,
…
,
T and save them; the last iteration of the filter gives us the probability distribution f
(
S
T
∣
ψ
T
), from which S
T is generated.Step 2.
Note that(40)
f
(
S
t
∣
S
t
+
1
,
ψ
t
)
=
f
(
S
t
,
S
t
+
1
∣
ψ
t
)
f
(
S
t
+
1
∣
ψ
t
)
=
f
(
S
t
+
1
∣
S
t
,
ψ
t
)
f
(
S
t
∣
ψ
t
)
f
(
S
t
+
1
∣
ψ
t
)
=
f
(
S
t
+
1
∣
S
t
)
f
(
S
t
∣
ψ
t
)
f
(
S
t
+
1
∣
ψ
t
)
∝
f
(
S
t
+
1
∣
S
t
)
f
(
S
t
∣
ψ
t
)
,
where f
(
S
t
+
1
∣
S
t
) is the transition probability and f
(
S
t
∣
ψ
t
) has been saved from Step 1. So we can generate S
t in the following way, first calculate
(41)
P
(
S
t
=
1
∣
S
t
+
1
,
ψ
t
)
=
f
(
S
t
+
1
∣
S
t
=
1
)
f
(
S
t
=
1
∣
ψ
t
)
∑
j
=
0
1
‍
f
(
S
t
+
1
∣
S
t
=
j
)
f
(
S
t
=
j
∣
ψ
t
)
,
and then generate S
t using a uniform distribution. For example, we generate a random number from a uniform distribution between 0 and 1; if this number is less than or equal to the calculated value of P
(
S
t
=
1
∣
S
t
+
1
,
ψ
t
), we set S
t
=
1, otherwise, S
t is set equal to 0.In view of applications, let us now consider the followingfour state MSM:
(42)
y
t
~
N
0
,
σ
S
t
2
,
t
=
1,2
,
…
,
T
,
σ
S
t
=
σ
1
S
1
,
t
+
σ
2
S
2
,
t
+
σ
3
S
3
,
t
+
σ
4
S
4
,
t
,
l
S
t
∈
{
1,2
,
3,4
}
with
transition
probabilities
p
i
j
:
=
P
(
S
t
=
j
∣
S
t
-
1
=
i
)
,
i
,
j
=
1,2
,
3,4
,
where S
k
,
t
=
1 if S
t
=
k, otherwise S
k
,
t
=
0. Note that this is a particular case of the model analysed in Section 2.1, where μ
t
=
0  ∀
t, hence we can perform the procedure referred to serially uncorrelated data taking μ
S
t
=
μ
=
0 to start the Gibbs sampling algorithm, therefore we haveStep 1.
GenerateS
~
T
:
=
(
S
1
,
S
2
,
…
,
S
T
) conditional on σ
~
2
:
=
(
σ
1
2
,
σ
2
2
,
σ
3
2
,
σ
4
2
),
(43)
p
~
:
=
p
11
,
p
12
,
p
13
,
p
21
,
p
22
,
p
23
,
p
31
,
p
32
,
p
33
,
p
41
,
p
42
,
p
43
,
ψ
T
:
=
y
1
,
y
2
,
…
,
y
T
.
For this purpose, we employ the Multi-move Gibbs sampling algorithm:(1)
run procedure in Section2.1 with μ
S
t
=
0 in order to get, from last iteration, f
(
S
T
∣
ψ
T
)
=
P
(
S
T
∣
ψ
T
),
(2)
recalling thatf
(
S
t
∣
S
t
+
1
,
ψ
t
)
∝
f
(
S
t
+
1
∣
S
t
)
f
(
S
t
∣
ψ
t
) for t
=
T
-
1
,
…
,
1, we can generate S
t from the vector of probabilities
(44)
P
(
S
t
=
1
∣
S
t
+
1
,
ψ
t
)
,
P
(
S
t
=
2
∣
S
t
+
1
,
ψ
t
)
,
h
P
(
S
t
=
3
∣
S
t
+
1
,
ψ
t
)
,
P
(
S
t
=
4
∣
S
t
+
1
,
ψ
t
)
,where, for i
=
1
,
…
,
4,
(45)
P
S
t
=
i
S
t
+
1
,
ψ
t
=
f
S
t
+
1
S
t
=
1
f
S
t
=
i
ψ
t
∑
j
=
1
3
‍
f
S
t
+
1
S
t
=
j
f
S
t
=
j
ψ
t
.Step 2.
Generateσ
~
2, conditional on S
~
T and the data ψ
T.
We want to impose the constraintσ
1
2
<
σ
2
2
<
σ
3
2
<
σ
4
2, so we redefine σ
S
t
2 in this way:
(46)
σ
S
t
2
=
σ
1
2
1
+
S
2
,
t
h
2
1
+
S
3
,
t
h
2
1
+
S
3
,
t
h
3
×
(
1
+
S
4
,
t
h
2
)
(
1
+
S
4
,
t
h
3
)
(
1
+
S
4
,
t
h
4
)
,
where h
j
>
0 for j
=
1
,
…
,
4, so that σ
2
2
:
=
σ
1
2
(
1
+
h
2
), σ
3
2
:
=
σ
1
2
(
1
+
h
2
)
(
1
+
h
3
) and σ
4
2
:
=
σ
1
2
(
1
+
h
2
)
(
1
+
h
3
)
(
1
+
h
4
). With this specification, we first generate σ
1
2, then generate h
-
2
=
1
+
h
2, h
-
3
=
1
+
h
3 and h
-
4
=
1
+
h
4 to obtain σ
2
2, σ
3
2 and σ
4
2 indirectly.Generating
σ
1
2
, Conditional on
h
2
,
h
3
and
h
4. Define for t
=
1
,
…
,
T
(47)
Y
t
1
:
=
y
t
(
1
+
S
2
,
t
h
2
)
(
1
+
S
3
,
t
h
2
)
(
1
+
S
3
,
t
h
3
)
(
1
+
S
4
,
t
h
2
)
(
1
+
S
4
,
t
h
3
)
(
1
+
S
4
,
t
h
4
)and takeY
1
t
~
N
(
0
,
σ
1
2
), in (42). By choosing aninverted Gamma prior distribution, that is, f
(
σ
1
2
∣
h
2
,
h
3
,
h
4
)
~
IG
(
ν
1
/
2
,
δ
1
/
2
) where ν
1
,
δ
1 are the known prior hyperparameters, it can be shown that the conditional posterior distribution from which we generate σ
1
2 is given by:
(48)
σ
1
2
∣
ψ
T
,
S
~
T
,
h
2
,
h
3
,
h
4
~
IG
ν
1
+
T
2
,
δ
1
+
∑
t
=
1
T
‍
Y
t
1
2
.Generating
h
-
2
Conditional on
σ
1
2
,
h
3
and
h
4. Note that the likelihood function of h
2 depends only on the values of y
t for which S
t
∈
{
2,3
,
4
}. Therefore, take y
t
(
1
)
:
=
{
y
t
∣
S
t
∈
{
2,3
,
4
}
,
t
=
1
,
…
,
T
} and denote with T
2 the size of this sample. Then define
(49)
Y
t
2
:
=
y
t
(
1
)
σ
1
2
(
1
+
S
3
,
t
h
3
)
(
1
+
S
4
,
t
h
3
)
(
1
+
S
4
,
t
h
4
)
,
hence, for the observation in which S
t
=
2, 3 or 4, we have Y
2
t
~
N
(
0
,
h
-
2
). If we choose an inverted Gamma distribution with parameters ν
2
,
δ
2 for the prior, we obtain h
-
2
=
1
+
h
2 from the following posterior distribution:
(50)
h
-
2
∣
ψ
T
,
S
~
T
,
σ
1
2
,
h
3
,
h
4
~
IG
ν
2
+
T
2
2
,
δ
2
+
∑
t
=
1
T
2
‍
Y
t
2
2
.
In case h
-
2
>
1 put h
2
=
h
-
2
-
1 and σ
2
2
:
=
σ
1
2
(
1
+
h
2
). Otherwise reiterate this step.Generating
h
-
3
Conditional on
σ
1
2
,
h
2
and
h
4. Operate in a similar way as above. In particular if we define y
t
(
2
)
:
=
{
y
t
∣
S
t
∈
{
3,4
}
,
t
=
1
,
…
,
T
} we will obtain
(51)
Y
t
3
:
=
y
t
(
2
)
σ
1
2
(
1
+
S
3
,
t
h
2
)
(
1
+
S
4
,
t
h
2
)
(
1
+
S
4
,
t
h
4
)
~
N
(
0
,
h
-
3
)
.Generating
h
-
4
Conditional on
σ
1
2
,
h
2
and
h
3. Operate in a similar way as above. In particular if we define y
t
(
3
)
:
=
{
y
t
∣
S
t
=
4
,
t
=
1
,
…
,
T
} we will have
(52)
Y
t
4
:
=
y
t
(
3
)
σ
1
2
(
1
+
S
4
,
t
h
2
)
(
1
+
S
4
,
t
h
3
)
~
N
(
0
,
h
-
4
)
.Step 3.
Generatep
~ conditional on S
~
T. In order to generate the transition probabilities we exploit the properties of theprior Beta distribution. Let us first define
(53)
p
-
i
i
:
=
P
(
S
t
≠
i
∣
S
t
-
1
=
i
)
=
1
-
p
i
i
,
i
=
1,2
,
3,4
;
p
-
i
j
:
=
P
(
S
t
=
j
∣
S
t
-
1
=
i
,
S
t
≠
i
)
,
i
≠
j
.
Hence we have that
(54)
p
i
j
:
=
P
(
S
t
=
j
∣
S
t
-
1
=
i
)
=
P
(
S
t
=
j
∣
S
t
-
1
=
i
,
S
t
≠
i
)
P
(
S
t
≠
i
∣
S
t
-
1
=
i
)
=
p
-
i
j
(
1
-
p
i
i
)
∀
i
≠
j
.
Given S
~
T, let n
i
j, i
,
j
=
1,2
,
3,4 be the total number of transitions from state S
t
-
1
=
i to S
t
=
j,  t
=
2,3
,
…
,
T and n
-
i
j the number of transitions from state S
t
-
1
=
i to S
t
≠
j.
Begin with the generation of probabilitiesp
i
i,  i
=
1,2
,
3,4 by taking the Beta distribution as conjugate prior: if we take p
i
i
~
Beta
(
u
i
i
,
u
-
i
i
), where u
i
i and u
-
i
i are the known hyperparameters of the priors, the posterior distribution of p
i
i given S
~
T still belongs to the Beta family distributions, that is,
(55)
p
i
i
∣
S
~
T
~
Beta
(
u
i
i
+
n
i
i
,
u
-
i
i
+
n
-
i
i
)
,
i
=
1,2
,
3,4
.
The others parameters, that is, p
i
j for i
≠
j and j
=
1,2
,
3, can be computed from the above equation p
i
j
=
p
-
i
j
(
1
-
p
i
i
), where p
-
i
j are generated from the following posterior Beta distribution:
(56)
p
-
i
j
∣
S
~
T
~
Beta
(
u
i
j
+
n
i
j
,
u
-
i
j
+
n
-
i
j
)
.
For example, given that p
11 is generated, we can obtain p
12 and p
13 by considering
(57)
p
-
12
∣
S
~
T
~
Beta
(
u
12
+
n
12
,
u
-
12
+
n
-
12
)
,
p
-
13
∣
S
~
T
~
Beta
(
u
13
+
n
13
,
u
-
13
+
n
-
13
)
,
where n
-
12
=
n
13
+
n
14 and n
-
13
=
n
12
+
n
14. Finally, given p
11, p
12 and p
13 generated in this way, we have p
14
=
1
-
p
11
-
p
12
-
p
13.Remark 2.
When we do not have any information about priors distribution we employ hyperparametersu
i
j
=
0.5, i
,
j
=
1,2
,
3,4. Usually we know that elements of the matrix diagonal in the transition matrix are bigger than elements out of the diagonal, because in a financial framework regime switching happens only occasionally: in this case, since we want p
i
i close to 1 and p
i
j, i
≠
j, close to 0, we will choose u
i
i bigger than u
-
i
i.


## 4. Goodness of Fit

Since financial time series are characterized by complex and rather unpredictable behavior, it is difficult to find, if there is any, a possible pattern. A typical set of techniques which allow tomeasure the goodness of forecasts obtained by using a certain model, is given by theresidual analysis. Let us suppose that we are provided with a time series of return observations {
y
t
}
t
=
1
,
…
,
T, T
∈
N
+, for which we choose, for example, the model described in (4) with M
=
4. By running the procedure of Section 2.1 we obtain the filtered probabilities
(58)
P
(
S
t
=
j
∣
ψ
t
)
,
j
=
1,2
,
3,4
,
t
=
1
,
…
,
T
,
and, by maximization of the log-likelihood function, we compute the parameters μ
^
1
,
…
,
μ
^
4
,
σ
^
1
,
…
,
σ
^
4, therefore we can estimate both the mean and variance of the process at time t, for any t
=
1
,
…
,
T, given the information set ψ
t as weighted average of four values:
(59)
μ
^
t
:
=
E
μ
t
∣
ψ
t
:
=
μ
^
1
P
S
t
=
1
∣
ψ
t
+
⋯
+
μ
^
4
P
(
S
t
=
4
∣
ψ
t
)
,
σ
^
t
2
:
=
E
(
σ
t
2
∣
ψ
t
)
:
=
σ
^
1
2
P
(
S
t
=
1
∣
ψ
t
)
+
⋯
+
σ
^
4
2
P
(
S
t
=
4
∣
ψ
t
)
.
If the chosen model fits well the data, then the standardized residuals will have the following form:
(60)
ϵ
^
t
=
y
t
-
μ
^
t
σ
^
t
~
N
0,1
,
t
=
1
,
…
,
T
,
therefore it is natural to apply a normality test as, for example, theJarque-Bera test (see [18]) for details. We recall briefly that Jarque-Bera statistics is defined as
(61)
JB
:
=
T
6
S
2
-
1
4
K
-
3
2
,
where the parameters S and K indicate the skewness, respectively, the kurtosis of ϵ
^
t. If ϵ
^
t come from a Normal distribution, the Jarque-Bera statistics converges asymptotically to a chi-squared distribution with two degrees of freedom, and can be used to test the null hypothesis of normality: this is a joint hypothesis of the skewness being zero and the excess kurtosis (
K
-
3
) being also zero.Remark 3.
Note that the Jarque-Bera test is very sensitive and often rejects the null hypothesis only because of a few abnormal observations, this is the reason why, one has to take point out these outliers which has to be canceled out before apply the test on the obtained smoothed data.


## 5. Prediction

The forecasting task is the most difficult step in the whole MSM approach. Let us suppose that our time series ends at timeT
∈
N
+, without further observations, then we have to start the prediction with the following quantities:(i)
the transition probability matrixP
*
:
=
{
p
i
j
}
i
,
j
=
1,2
,
3,4;
(ii)
the vectorπ
T
:
=
P
(
S
T
∣
ψ
T
)
=
(
P
(
S
T
=
1
∣
ψ
T
)
,
…
,
P
(
S
T
=
4
∣
ψ
T
)
) obtained from the last iteration of the filter algorithm, for example, the procedure in Section 2.1.It follows that we have to proceed with the first step of the filter procedure obtaining the one-step ahead probability of the state S
T
+
1 given the sample of observations ψ
T, that is,
(62)
P
(
S
T
+
1
=
j
∣
ψ
T
)
=
∑
i
=
1
4
‍
p
i
j
P
(
S
T
=
j
∣
ψ
T
)
,
j
=
1,2
,
3,4
.
Equation (62) can be seen as a prediction for the regime at time T
+
1, knowing observations up to time T. At this point, the best way to make prediction about the unobserved variable is the simulation of further observations. Indeed, with the new probability P
(
S
T
+
1
∣
ψ
T
) and the vector of parameter estimates θ
^
:
=
(
μ
^
1
,
…
,
μ
^
4
,
σ
^
1
,
…
,
σ
^
4
) we can estimate the one step ahead mean and variance as follows:
(63)
μ
^
T
+
1
:
=
E
(
μ
T
+
1
∣
ψ
T
)
:
=
μ
^
1
P
(
S
T
+
1
=
1
∣
ψ
T
)
+
⋯
+
μ
^
4
P
(
S
T
+
1
=
4
∣
ψ
T
)
,
σ
^
T
+
1
2
:
=
E
(
σ
T
+
1
2
∣
ψ
T
)
:
=
σ
^
1
2
P
(
S
T
+
1
=
1
∣
ψ
T
)
+
⋯
+
σ
^
4
2
P
(
S
T
+
1
=
4
∣
ψ
T
)
.
Then we simulate y
^
T
+
1 by the Gaussian distribution N
(
μ
^
T
+
1
,
σ
^
T
+
1
) and, once y
^
T
+
1 has been simulated, we defe ψ
T
+
1
:
=
{
y
1
,
…
,
y
T
,
y
^
T
+
1
}. Then we first apply again the filter procedure of Section 2.1 for t
=
T
+
1 in order to obtain P
(
S
T
+
1
∣
ψ
T
+
1
), then we compute P
(
S
T
+
1
∣
ψ
T
+
1
),  μ
^
T
+
2 and σ
^
T
+
2
2, and we simulate y
^
T
+
2 by the Gaussian distribution N
(
μ
^
T
+
2
,
σ
^
T
+
2
). Latter procedure runs the same all the other rime-steps T
+
3
,
…
,
T
+
m, where m
∈
N
+ is the time horizon of our forecast.Remark 4.
We would like to underline that latter described method is not reliable with few simulations since eachy
T
+
τ, for τ
=
1
,
…
,
m may assume a wide range of values and a single drawn describes only one of the many possible paths. So we can think to reiterate previous strategy many times in order to compute the mean behavior of P
(
S
T
+
τ
∣
ψ
T
+
τ
), μ
^
T
+
τ and σ
^
T
+
τ. After having obtained a satisfactory number of data, then we can construct a confidence interval within the state probability will more likely take value. Obviously a high number of iterations of latter procedure rapidly increases the computational complexity of the whole algorithm because of the MLE related computational complexity, therefore we will adopt a rather different strategy which consists in simulating y
T
+
τ
N times at each step (e.g., N
=
10000) and then taking the mean over those values. However, we must pay attention because the mean calculation could cancel the possible regime switching: for example, if we draw many times y
t from N
(
0
,
σ
S
t
) and we take the mean, by the law of large number we will have zero at any time. To overcome this problem we can take the mean of absolute values and then multiply this mean by a number x, which is a random variable that takes values 1 or −1, with equal probability, hence deciding the sign of y
t at every simulation step.


## 6. Applications

In this section we are going to apply the classical inference approach for a MSM to analyse real financial time series. In particular we will first examine data coming from the Standard & Poor’s 500 (S&P 500) equity index which is considered, being based on the 500 most important companies in the United States, as one of the best representations of the U.S. stock market. Secondly, we shall consider the DAX (Deutsche Aktien Index) index which follows the quotations of the 30 major companies in Germany. Our choice is motivated by a twofold goal: first we want to test the proposed 4-states MSM model on two particularly significant indexes which have shown to incorporate abrupt changes and oscillations, secondly we aim at comparing the behaviour of the two indexes between each other.Computations have been performed following the MSM approach described in previous section, namely exploiting the procedures illustrated in Section2. Let us underline that, instead of a standard 3-states MSM model, we shall use a 4-states MSM approach both for the S&P 500 and the DAX returns. Moreover the analysis has been realized for different intervals of time, focusing mainly on the period of Global Financial Crisis.


### 6.1. The S&P 500 Case

Figure1 reports the graph of the Standard & Poor’s 500 from 1st June 1994 to 27th May 2014, and it points out the dramatic collapse of index prices in years 2008-2009, when the crisis blowed-up causing the achievement, 6th of March 2009 with 683.38 points, of the lowest value since September 1996.Figure 1
Daily observations of S&P 500 from 1994 to 2014.Because of the latter fact we decided to focus our analysis on recent years. In particular we take into account data starting from the 1st of June 2007, and until 27 May 2014, therefore, denoting withΛ the set of observations and with X
t,  t
∈
Λ, the price data of the S&P 500, returns are calculated as usual by y
t
:
=
(
X
t
-
X
t
-
1
)
/
X
t
-
1, t
∈
Λ, where {
y
t
}
t
∈
Λ are the values for which we want to choose thebest MSM. Note that in our implementation we grouped the daily data into weekly parcels in order to make the filter procedures less time-consuming and have a more clear output, therefore we obtain a vector of 365 values, still denoted by y
t, as shown in Figure 2.Figure 2
Daily returns of S&P 500 from 2007 to 2014.Next step consists in understand if the returns are serially correlated or serially uncorrellated, a taks which can be accomplished by running some suitable test, for example, the Durbin-Watson test (see, for example, [19, 20] or [7]) computing directly the value of the autoregressive parameter ϕ by least square methods, namely ϕ
:
=
(
∑
t
∈
Λ
‍
y
t
y
t
+
1
)
/
(
∑
t
∈
Λ
‍
y
t
2
), which gives us a rather low value, that is, -
0.0697, so that we can neglect the autoregressive pattern and start the analysis by considering S&P 500 returns to be generated by a Gaussian distribution with switching mean and variance, that is,
(64)
y
t
=
μ
S
t
+
ϵ
t
,
t
∈
Λ
ϵ
t
t
∈
{
1
,
…
,
T
}
are
i
.
i
.
d
.
N
(
0
,
σ
S
t
2
)
,
μ
S
t
=
μ
1
S
1
,
t
+
⋯
+
μ
4
S
4
,
t
,
σ
S
t
=
σ
1
S
1
,
t
+
⋯
+
σ
4
S
4
,
t
,
S
t
∈
{
1
,
…
,
4
}
with
probabilities
p
i
j
,
where for (
k
,
t
)
∈
{
1
,
…
,
4
}
×
Λ we have S
k
,
t
=
1 if S
t
=
k, otherwise S
k
,
t
=
0. Therefore, we suppose that the state variable S
t, t
∈
Λ, takes its values in the set Ω
:
=
{
1,2
,
3,4
}, and we expect that the probabilities of being in the third and fourth state increase as a financial crisis occurs. Exploiting the procedure provided in Section 2.1, with respect to the returns y
t, t
∈
Λ
=
{
1
,
…
,
365
}, we get the results shown in Figures 3 and 4.Figure 3
Filtered probabilitiesP
(
S
t
=
k
∣
ψ
t
) for t
∈
Λ, k
∈
{
1,2
,
3,4
}.Figure 4
Smoothed probabilitiesP
(
S
t
=
k
∣
ψ
T
) for t
∈
Λ, k
∈
{
1,2
,
3,4
}.Let us now consider the estimated standard deviation(65)
σ
^
t
:
=
E
(
σ
t
∣
ψ
t
)
:
=
σ
^
1
P
(
S
t
=
1
∣
ψ
t
)
+
⋯
+
σ
^
4
P
(
S
t
=
4
∣
ψ
t
)
,
t
∈
Λ
,
which we want to compare with the VIX index, also known as theChicago Board Options Exchange (CBOE) market volatility index, namely one of the most relevant measure for the implied volatility of the S&P 500 index, whose value used by our analysis are reported in Figure 5.Figure 5
CBOE volatility index (VIX), daily data from 1994 to 2014.What we obtain by plotting both estimated volatility and VIX values in the same graph can be seen in Figure6, where the VIX trend is plotted in red, while we have used the blue color for the conditional standard deviation values.Figure 6
VIX index (red) versus estimated volatility (blue).Note that, in order to have values of the same order, eachσ
^
t, t
∈
Λ, has been multiplied by a scale factor equal to 1000. We would like to point out how the estimated standard deviation accurately approximates the VIX behaviour, hence allow us to count on an effectivesubstitute for the volatility of the S&P 500, at least during a relative nonchaotic period. In fact we also underline that the greater discrepancies between real and simulated values appears during the maximum intensity period of the recent financial crisis. In particular the widest gaps are realized in correspondence with the recession experienced at the end of 2008.In what follows we study how latter evidence influences the global goodness of the realized analysis. In particular we performed agoodness of fit analysis computing the standardized residuals of the proposed MSM by ϵ
^
t
:
=
(
y
t
-
μ
^
t
)
/
σ
^
t, t
∈
Λ, where y
t is the observation of S&P 500 return at time t, μ
^
t is the estimated conditional mean, and σ
^
t is the standard deviation. If the model is a good fit for the S&P 500 return, standardized residuals will be generated by a standard Gaussian distribution. In Figures 7 and 8, we have reported both the histogram, its related graph, and the so callednormal probability plot (NPP) for the standardized residuals.Figure 7
Plot and histogram of standardized residuals.Figure 8
Normal probability plot of standardized residuals.Let us recall that the purpose of the NPP is to graphically assess whether the residuals could come from a normal distribution. Indeed, if such a hypothesis holds, then theNPP has to be linear, namely the large majority of the computed values, that is, the blue points in Figure 8, should stay close to aparticular line, which is the red dotted one in Figure 8, which is the case in our analysis apart from the three points in the left-hand corner of the graph, which correspond to the minimal values of the vector of standardized residuals.Applying two normality tests on{
ϵ
^
t
}
t
∈
Λ, that is, theJarque-Bera test and (see, for example, [21, pag.  443]) theLilliefors test, we have that the null hypothesis of normality for the standardized residuals can be rejected at the 5% level, unless the previous pointed outoutliers are removed. Indeed if the two minimal standardized residuals, corresponding to ϵ
^
71
=
-
3.8441 and ϵ
^
153
=
-
3.6469, are cancelled out from the vector {
ϵ
^
t
}
t
∈
Λ, previously cited tests indicate that the hypothesis of normality, at the same significance level of 5
%, cannot be rejected. In particular, the Jarque-Bera statistics value is JB
=
2.7858 with corresponding P-value P
JB
=
0.2153, and the critical value for this test, that is, the maximum value of the JB statistics for which the null hypothesis cannot be rejected at the chosen significance level, is equal to k
JB
=
5.8085. Similarly, with regard to the Lilliefors test, numerical value of Liellifors statistics, P-value and critical value are respectively given by L
=
0.0424, P
L
=
0.1181 and k
L
=
0.0472.In what follows we develop the forecast procedure shown in Section5. Since we are dealing with weekly data, let us suppose we want to predict probability of volatility σ
^
t, t
∈
Λ on a time horizon of two months, hence 8 steps ahead, then simulations have been performed according to Remark 4, with N
=
15000, T
=
365, τ
=
{
1,2
,
…
,
8
}, and x uniformly distributed in {
-
1,1
}. Obtained forecasting results are shown in Figure 9, where plots are referred to the observations from the 300th to the 373rd, with the last 8 simulated values within red rectangles.Forecast for returns (a) and conditional standard deviation (b).
(a)
(b)


### 6.2. The DAX Case

In what follows the proposed 4-state MSM shall be applied to analyse theDeutsche Aktien Index (DAX) stock index during a shorter, compared to the study made in Section 6.1, time interval, indeed we take into account data between the 1st of January 2008 and until the 1st of January 2011. Such a choice is motivated by the will of concentrate our attention on the most unstable period, that for the DAX index starts with the abrupt fall experienced by the European markets at the beginning of 2008. We underline that, with respect to the considered time period, the estimated autoregressive coefficient for the DAX index is a bit higher, in absolute value, than the one obtained for the S&P 500 index, indeed ϕ
^
DAX
=
-
0.1477 versus ϕ
^
S
&
P
500
=
-
0.0697, so that we decided to fit the weekly DAX returns with an autoregressive pattern using the procedure in Section 2.2 and according to the following MSM:
(66)
y
t
=
ϕ
(
y
t
-
1
)
+
ϵ
t
,
t
∈
Λ
ϵ
t
t
∈
{
1
,
…
,
T
}
are
i
.
i
.
d
.
N
(
0
,
σ
S
t
2
)
,
σ
S
t
=
σ
1
S
1
,
t
+
⋯
+
σ
4
S
4
,
t
,
S
t
∈
{
1
,
…
,
4
}
with
probabilities
p
i
j
,
where Λ
=
{
1,2
,
…
,
156
} is the set of times, y
t is the return of DAX index at time t and for (
k
,
t
)
∈
{
1
,
…
,
4
}
×
Λ we have S
k
,
t
=
1 if S
t
=
k, otherwise S
k
,
t
=
0. Note that we have slightly simplified the model keeping μ
S
t
=
0, that is we are dealing with a switching variance model where only σ
S
t switches in regime. The latter is not a restriction since we are interested in model the DAX index volatility. In Figure 10 we report the plot of weekly returns within the period of interest.Figure 10
Weekly returns of DAX from 1st January 2008 to 1st January 2011.Figures11, 12 and 13, show the graphical outputs of the performed analysis. In particular they respectively display both the filtered and smoothed probabilities for every regime, and estimated standard deviation σ
^
t, for t
∈
Λ.Figure 11
Filtered probabilitiesP
(
S
t
=
k
∣
ψ
t
) for t
∈
Λ, k
∈
{
1,2
,
3,4
}.Figure 12
Smoothed probabilitiesP
(
S
t
=
k
∣
ψ
T
) for t
∈
Λ, k
∈
{
1,2
,
3,4
}.Figure 13
Estimated conditional standard deviation for DAX index.Thegoodness of fit analysis has been done applying normality tests to standardized residuals, with related NPP shown in Figure 14.Figure 14
Normal probability plot for standardized residuals.As in the S&P 500 case (see Section6.1) both the Jarque-Bera and Lilliefors tests reject the null hypothesis of normality if they are performed taking into account the whole set of residuals, while, removing the single outlier corresponding to the smallest value of residuals (see the left-hand corner of the NPP in Figure 14) both tests are passed and we can conclude that our 4-state MSM provides a good fitting of examined data. In particular, as in the analysis of S&P 500, we computed numerical values for the two statistical test with a significance level α
=
5
%, obtaining JB
=
0.2326, P
JB
=
0.500, k
JB
=
5.6006 and L
=
0.0712, P
L
=
0.0553, k
L
=
0.0720.


## 6.1. The S&P 500 Case

Figure1 reports the graph of the Standard & Poor’s 500 from 1st June 1994 to 27th May 2014, and it points out the dramatic collapse of index prices in years 2008-2009, when the crisis blowed-up causing the achievement, 6th of March 2009 with 683.38 points, of the lowest value since September 1996.Figure 1
Daily observations of S&P 500 from 1994 to 2014.Because of the latter fact we decided to focus our analysis on recent years. In particular we take into account data starting from the 1st of June 2007, and until 27 May 2014, therefore, denoting withΛ the set of observations and with X
t,  t
∈
Λ, the price data of the S&P 500, returns are calculated as usual by y
t
:
=
(
X
t
-
X
t
-
1
)
/
X
t
-
1, t
∈
Λ, where {
y
t
}
t
∈
Λ are the values for which we want to choose thebest MSM. Note that in our implementation we grouped the daily data into weekly parcels in order to make the filter procedures less time-consuming and have a more clear output, therefore we obtain a vector of 365 values, still denoted by y
t, as shown in Figure 2.Figure 2
Daily returns of S&P 500 from 2007 to 2014.Next step consists in understand if the returns are serially correlated or serially uncorrellated, a taks which can be accomplished by running some suitable test, for example, the Durbin-Watson test (see, for example, [19, 20] or [7]) computing directly the value of the autoregressive parameter ϕ by least square methods, namely ϕ
:
=
(
∑
t
∈
Λ
‍
y
t
y
t
+
1
)
/
(
∑
t
∈
Λ
‍
y
t
2
), which gives us a rather low value, that is, -
0.0697, so that we can neglect the autoregressive pattern and start the analysis by considering S&P 500 returns to be generated by a Gaussian distribution with switching mean and variance, that is,
(64)
y
t
=
μ
S
t
+
ϵ
t
,
t
∈
Λ
ϵ
t
t
∈
{
1
,
…
,
T
}
are
i
.
i
.
d
.
N
(
0
,
σ
S
t
2
)
,
μ
S
t
=
μ
1
S
1
,
t
+
⋯
+
μ
4
S
4
,
t
,
σ
S
t
=
σ
1
S
1
,
t
+
⋯
+
σ
4
S
4
,
t
,
S
t
∈
{
1
,
…
,
4
}
with
probabilities
p
i
j
,
where for (
k
,
t
)
∈
{
1
,
…
,
4
}
×
Λ we have S
k
,
t
=
1 if S
t
=
k, otherwise S
k
,
t
=
0. Therefore, we suppose that the state variable S
t, t
∈
Λ, takes its values in the set Ω
:
=
{
1,2
,
3,4
}, and we expect that the probabilities of being in the third and fourth state increase as a financial crisis occurs. Exploiting the procedure provided in Section 2.1, with respect to the returns y
t, t
∈
Λ
=
{
1
,
…
,
365
}, we get the results shown in Figures 3 and 4.Figure 3
Filtered probabilitiesP
(
S
t
=
k
∣
ψ
t
) for t
∈
Λ, k
∈
{
1,2
,
3,4
}.Figure 4
Smoothed probabilitiesP
(
S
t
=
k
∣
ψ
T
) for t
∈
Λ, k
∈
{
1,2
,
3,4
}.Let us now consider the estimated standard deviation(65)
σ
^
t
:
=
E
(
σ
t
∣
ψ
t
)
:
=
σ
^
1
P
(
S
t
=
1
∣
ψ
t
)
+
⋯
+
σ
^
4
P
(
S
t
=
4
∣
ψ
t
)
,
t
∈
Λ
,
which we want to compare with the VIX index, also known as theChicago Board Options Exchange (CBOE) market volatility index, namely one of the most relevant measure for the implied volatility of the S&P 500 index, whose value used by our analysis are reported in Figure 5.Figure 5
CBOE volatility index (VIX), daily data from 1994 to 2014.What we obtain by plotting both estimated volatility and VIX values in the same graph can be seen in Figure6, where the VIX trend is plotted in red, while we have used the blue color for the conditional standard deviation values.Figure 6
VIX index (red) versus estimated volatility (blue).Note that, in order to have values of the same order, eachσ
^
t, t
∈
Λ, has been multiplied by a scale factor equal to 1000. We would like to point out how the estimated standard deviation accurately approximates the VIX behaviour, hence allow us to count on an effectivesubstitute for the volatility of the S&P 500, at least during a relative nonchaotic period. In fact we also underline that the greater discrepancies between real and simulated values appears during the maximum intensity period of the recent financial crisis. In particular the widest gaps are realized in correspondence with the recession experienced at the end of 2008.In what follows we study how latter evidence influences the global goodness of the realized analysis. In particular we performed agoodness of fit analysis computing the standardized residuals of the proposed MSM by ϵ
^
t
:
=
(
y
t
-
μ
^
t
)
/
σ
^
t, t
∈
Λ, where y
t is the observation of S&P 500 return at time t, μ
^
t is the estimated conditional mean, and σ
^
t is the standard deviation. If the model is a good fit for the S&P 500 return, standardized residuals will be generated by a standard Gaussian distribution. In Figures 7 and 8, we have reported both the histogram, its related graph, and the so callednormal probability plot (NPP) for the standardized residuals.Figure 7
Plot and histogram of standardized residuals.Figure 8
Normal probability plot of standardized residuals.Let us recall that the purpose of the NPP is to graphically assess whether the residuals could come from a normal distribution. Indeed, if such a hypothesis holds, then theNPP has to be linear, namely the large majority of the computed values, that is, the blue points in Figure 8, should stay close to aparticular line, which is the red dotted one in Figure 8, which is the case in our analysis apart from the three points in the left-hand corner of the graph, which correspond to the minimal values of the vector of standardized residuals.Applying two normality tests on{
ϵ
^
t
}
t
∈
Λ, that is, theJarque-Bera test and (see, for example, [21, pag.  443]) theLilliefors test, we have that the null hypothesis of normality for the standardized residuals can be rejected at the 5% level, unless the previous pointed outoutliers are removed. Indeed if the two minimal standardized residuals, corresponding to ϵ
^
71
=
-
3.8441 and ϵ
^
153
=
-
3.6469, are cancelled out from the vector {
ϵ
^
t
}
t
∈
Λ, previously cited tests indicate that the hypothesis of normality, at the same significance level of 5
%, cannot be rejected. In particular, the Jarque-Bera statistics value is JB
=
2.7858 with corresponding P-value P
JB
=
0.2153, and the critical value for this test, that is, the maximum value of the JB statistics for which the null hypothesis cannot be rejected at the chosen significance level, is equal to k
JB
=
5.8085. Similarly, with regard to the Lilliefors test, numerical value of Liellifors statistics, P-value and critical value are respectively given by L
=
0.0424, P
L
=
0.1181 and k
L
=
0.0472.In what follows we develop the forecast procedure shown in Section5. Since we are dealing with weekly data, let us suppose we want to predict probability of volatility σ
^
t, t
∈
Λ on a time horizon of two months, hence 8 steps ahead, then simulations have been performed according to Remark 4, with N
=
15000, T
=
365, τ
=
{
1,2
,
…
,
8
}, and x uniformly distributed in {
-
1,1
}. Obtained forecasting results are shown in Figure 9, where plots are referred to the observations from the 300th to the 373rd, with the last 8 simulated values within red rectangles.Forecast for returns (a) and conditional standard deviation (b).
(a)
(b)


## 6.2. The DAX Case

In what follows the proposed 4-state MSM shall be applied to analyse theDeutsche Aktien Index (DAX) stock index during a shorter, compared to the study made in Section 6.1, time interval, indeed we take into account data between the 1st of January 2008 and until the 1st of January 2011. Such a choice is motivated by the will of concentrate our attention on the most unstable period, that for the DAX index starts with the abrupt fall experienced by the European markets at the beginning of 2008. We underline that, with respect to the considered time period, the estimated autoregressive coefficient for the DAX index is a bit higher, in absolute value, than the one obtained for the S&P 500 index, indeed ϕ
^
DAX
=
-
0.1477 versus ϕ
^
S
&
P
500
=
-
0.0697, so that we decided to fit the weekly DAX returns with an autoregressive pattern using the procedure in Section 2.2 and according to the following MSM:
(66)
y
t
=
ϕ
(
y
t
-
1
)
+
ϵ
t
,
t
∈
Λ
ϵ
t
t
∈
{
1
,
…
,
T
}
are
i
.
i
.
d
.
N
(
0
,
σ
S
t
2
)
,
σ
S
t
=
σ
1
S
1
,
t
+
⋯
+
σ
4
S
4
,
t
,
S
t
∈
{
1
,
…
,
4
}
with
probabilities
p
i
j
,
where Λ
=
{
1,2
,
…
,
156
} is the set of times, y
t is the return of DAX index at time t and for (
k
,
t
)
∈
{
1
,
…
,
4
}
×
Λ we have S
k
,
t
=
1 if S
t
=
k, otherwise S
k
,
t
=
0. Note that we have slightly simplified the model keeping μ
S
t
=
0, that is we are dealing with a switching variance model where only σ
S
t switches in regime. The latter is not a restriction since we are interested in model the DAX index volatility. In Figure 10 we report the plot of weekly returns within the period of interest.Figure 10
Weekly returns of DAX from 1st January 2008 to 1st January 2011.Figures11, 12 and 13, show the graphical outputs of the performed analysis. In particular they respectively display both the filtered and smoothed probabilities for every regime, and estimated standard deviation σ
^
t, for t
∈
Λ.Figure 11
Filtered probabilitiesP
(
S
t
=
k
∣
ψ
t
) for t
∈
Λ, k
∈
{
1,2
,
3,4
}.Figure 12
Smoothed probabilitiesP
(
S
t
=
k
∣
ψ
T
) for t
∈
Λ, k
∈
{
1,2
,
3,4
}.Figure 13
Estimated conditional standard deviation for DAX index.Thegoodness of fit analysis has been done applying normality tests to standardized residuals, with related NPP shown in Figure 14.Figure 14
Normal probability plot for standardized residuals.As in the S&P 500 case (see Section6.1) both the Jarque-Bera and Lilliefors tests reject the null hypothesis of normality if they are performed taking into account the whole set of residuals, while, removing the single outlier corresponding to the smallest value of residuals (see the left-hand corner of the NPP in Figure 14) both tests are passed and we can conclude that our 4-state MSM provides a good fitting of examined data. In particular, as in the analysis of S&P 500, we computed numerical values for the two statistical test with a significance level α
=
5
%, obtaining JB
=
0.2326, P
JB
=
0.500, k
JB
=
5.6006 and L
=
0.0712, P
L
=
0.0553, k
L
=
0.0720.


## 7. Conclusion

In this work we have shown how a four-state Regime Switching Model can be successfully exploited to study the volatility parameter which strongly characterizes concrete financial time series. In particular we have performed a deep analysis of data taken from theStandard & Poor’s 500 and theDeutsche Aktien Index series of returns focusing our attention on the most erratic periods for financial quantities influenced by the global financial crisis. We provided both accurate numerical results and related effective graphs, together with a detailed overview of the algorithms used to infer on parameters of the exploited switching models. Applications to the above mentioned markets (see Section 6) are completed by an accurate good of fitness study.



---

*Source: 753852-2014-12-18.xml*