Markov-Switching, Copula model estimation - Ca' Foscari [PDF]

“Markov-Switching copula models for dependence analysis in time series”. Supervisor: Mr. Casarin. Student: Victor Fa

30 downloads 20 Views 3MB Size

Recommend Stories


Annali di Ca' Foscari
Be like the sun for grace and mercy. Be like the night to cover others' faults. Be like running water

università ca' foscari venezia
Be who you needed when you were younger. Anonymous

Ca' Foscari University Course Guide
Respond to every call that excites your spirit. Rumi

Ca' Foscari School for International Education in Venice
Be like the sun for grace and mercy. Be like the night to cover others' faults. Be like running water

Ca' Foscari School for International Education in Venice
We may have all come on different ships, but we're in the same boat now. M.L.King

Università Ca' Foscari, Venezia CORSO DI ALTA FORMAZIONE IN COMUNICAZIONE
The happiest people don't have the best of everything, they just make the best of everything. Anony

Model Fitting and Error Estimation [PDF]
A mathematical model designed to fit experimental data so as to explicitly .... Linear Regression Analysis. E(a,b) = ..... www.nr.com/nronline_switcher.html. • Matlab online help www.mathworks.com/access/helpdesk/help/techdoc/. • References. Ande

Is a normal copula the right copula?
Make yourself a priority once in a while. It's not selfish. It's necessary. Anonymous

Levy copula
Life isn't about getting and having, it's about giving and being. Kevin Kruse

Estimation of Non-Simplified Vines Using Nonparametric Trivariate Copula Constructions
If you want to become full, let yourself be empty. Lao Tzu

Idea Transcript


“Markov-Switching copula models for dependence analysis in time series”

Supervisor: Mr. Casarin Student: Victor Fardel

Master Degree, Financial Econometrics

Ca’ foscari University June 2014

2

CONTENTS

Contents 1 Introduction

5

2 Methodology

6

2.1

MCMC simulation method . . . . . . . . . . . . . . . . . . . . . . .

6

2.2

Bayesian estimation . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

2.3

Simulation and estimation processes . . . . . . . . . . . . . . . . . . 10

2.4

2.3.1

Gibbs sampling methodology . . . . . . . . . . . . . . . . . 10

2.3.2

Metropolis-Hasting algorithm . . . . . . . . . . . . . . . . . 12

2.3.3

Prior distributions of the parameters . . . . . . . . . . . . . 13

2.3.4

Hamilton filter . . . . . . . . . . . . . . . . . . . . . . . . . 16

Copula . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.4.1

Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

2.4.2

The Gaussian copula model . . . . . . . . . . . . . . . . . . 18

3 Data

20

3.1

Industrial Production Index . . . . . . . . . . . . . . . . . . . . . . 20

3.2

The electricity demand in Australia . . . . . . . . . . . . . . . . . . 24

4 Results

29

4.1

Simulation process . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

4.2

On the Industrial Production Index (IPI) . . . . . . . . . . . . . . . 36

4.3

On the Australian electricity market . . . . . . . . . . . . . . . . . 40

5 Conclusion

46

6 Annexe

48

3

LIST OF FIGURES 6.1

Matlab code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 6.1.1

DGP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

6.1.2

Gibbs sampling . . . . . . . . . . . . . . . . . . . . . . . . . 51

6.1.3

Forward Filtering Backward Sampling(FFBS) . . . . . . . . 55

6.1.4

Hamilton filter . . . . . . . . . . . . . . . . . . . . . . . . . 56

6.1.5

Parameter estimation . . . . . . . . . . . . . . . . . . . . . . 56

List of Tables 1

Descriptive statistics IPI. . . . . . . . . . . . . . . . . . . . . . . . . 22

2

Descriptive statistics of the growth rate of the IPI. . . . . . . . . . . 24

3

Descriptive statistics of the demand of electricity series. . . . . . . . 26

4

Descriptive statistics of the growth rate of the demand of electricity series. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

5

Convergence results. . . . . . . . . . . . . . . . . . . . . . . . . . . 35

6

Convergence results. . . . . . . . . . . . . . . . . . . . . . . . . . . 39

7

Convergence results. . . . . . . . . . . . . . . . . . . . . . . . . . . 45

List of Figures 1

Industrial Production Index (IPI). . . . . . . . . . . . . . . . . . . . 21

2

IPI monthly, quarterly and annually growth rate. . . . . . . . . . . 23

3

Daily demand of electricity in New South Wales and in Victoria state. 25

4

Daily growth rate of the demand of electricity in New South Wales and in Victoria state. . . . . . . . . . . . . . . . . . . . . . . . . . . 27

4

LIST OF FIGURES 5

Histogram of the daily growth rate of the demand of electricity in New South Wales and in Victoria state. . . . . . . . . . . . . . . . . 29

6

Simulated data and simulated and forecast state. . . . . . . . . . . 30

7

State value beyond the Gibbs sampling. . . . . . . . . . . . . . . . . 31

8

State value beyond the Gibbs sampling. . . . . . . . . . . . . . . . . 32

9

Conditional mean value beyond the Gibbs sampling. . . . . . . . . . 33

10

Conditional variance value beyond the Gibbs sampling. . . . . . . . 34

11

Estimation of the correlation coefficient, copula part . . . . . . . . . 35

12

Growth rate of the IPI in the Euro Area and in the USA, forecast state and probabilities. . . . . . . . . . . . . . . . . . . . . . . . . . 36

13

Estimate state of the IPI series. . . . . . . . . . . . . . . . . . . . . 37

14

Estimation of the probabilities in each state for the two EA and the USA. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

15

Estimation of the mean in each state for the two EA and the USA.

38

16

Estimation of the variance in each state for the two EA and the USA. 39

17

Estimation of the correlation coefficient, copula part. . . . . . . . . 40

18

Growth rate of the electricity demand, forecast state and probabilities. 41

19

Estimation of the probabilities in each state for the NSW and the Victoria states. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

20

Estimation of the mean in each state for the two EA and the USA.

43

21

Estimation of the variance in each state for the two EA and the USA. 44

22

Estimation of the correlation coefficient between, copula part. . . . 45

5

1

1 INTRODUCTION

Introduction

In the thesis, we would like to study the behavior of the series of interest. First of all, we would like to be able to describe the state of the series (is it a period of expansion or a recession period, is there a high demand or a low demand for a specific commodity). After that, we would like to evaluate the correlation between the two series of interest. After applying our model over simulated data, we will go through observed data sets. On one hand, we will study the Industrial Production Index (IPI) for the Euro Area and the United States of America. On a other hand, we will study the Australian electricity market, especially the demand for electricity in the state of Victoria and Queensland. In both of the application case, the study is interesting: on the IPI market, the results will provide us with the dependence between the European and the American economies: are they correlated, and if they are what about the correlation? On the electricity market, the results will provide us also with the dependence between the two series. Our results on the IPI could be interesting on an political point of view. They provide the politician with the different states of the economy and the dependence between the two economic areas. For example, the politician could be interested to compare the characteristics of the economy in the different states. Plus, the measure of dependence will inform them about a dependence between the economic activity of the two economic area and provide them with a measure of dependence of the two economy. Based on the assumption made on our models, two states for the Markov-Switching process, the same prior distribution for variable of interest, we would expected to find convergent estimation of all our parameter: mean, variance, state probabilities and correlation coefficient. We would be able to provide a coherent description of the each state of the variables of interest and a measure of dependence between two variables. The remainder of the thesis will be organized as follow: in a first part, we will go

6

2 METHODOLOGY

through the methodology we will use in the thesis. Afterwards, in a second part, we will present the data samples (the simulated data sample and the observed data on the Industrial Production Index (IPI) and the electricity growth rate of the demand). To end, in a third part, we will provide you with the empirical results both on the state estimation for each variable and the correlation between the variables.

2

Methodology

As explained in the introduction, the aim of the thesis is to provide a method to point out the relation, or the stochastic dependence, between two time series. We illustrate the effectiveness of the proposed method with two examples. In the first case we will apply our methodology over the Industrial Production Index (IPI) to study the relationship between the level of the economic activity in the United States of America (USA) and the value for the Euro Area (EA). In the second example, we will go through a financial dataset and we will apply our methodology on the Australian electricity market data in the aim to observe the regional dependence of the series. Before that, we would like to be able to determine the regime in which the returns of an asset are.

2.1

MCMC simulation method

Our first objective is to provide a flexible model for univariate time series. Actually, we can express the returns of time series as given by the equation 2: a conditional mean (µ) and a error term (). The data meander around the mean.

rt = log(yt ) − log(yt−1 )

(1)

rt = µt + t

(2)

7

2 METHODOLOGY

Where yt are the raw data, and t is the time index. Note that the logarithmic difference provides us with an asymptotic estimator of the growth rate of the data. However, as we know, following the idea express in the article [Engle et Ng(1993)], positive and negative has not the same impact on the series. Most of the time, negative shocks would imply more volatility than positive shocks. Because of that, we have chosen to define the return series as Markov Switching (MS) model as define thanks to the equation 3. We have chosen to use arbitrary two different regimes, then we allowed the return to have different means and variance according to the state.

rt =

                        

µ1 + t , St = 1 t ∼ N (0, σ12 ) (3) µ2 + t , St = 2 t ∼ N (0, σ22 )

Note that, St is a latent state variable, an unobservable variable. In our case, this variable will take two values 1 and 2. According to the prior condition, regime 1 and regime 2 will correspond to refer to different economic situation. Thus, we can express the model as follow:

rt | St ∼ N (µ(St ), σ(St )) µ(St ) =

σ(St ) =

  

µ1 , if St = 1

    

µ2 , if St = 2

 

σ2 , if St = 2

σ1 , if St = 1

Plus, St in our case follows a Markov Chain (MC). A MC is defined as a process for which the present realization of the process only depends on the most recent

8

2 METHODOLOGY

past realization. The equation 4 provide us with the transition probability of a N-states MC. pij refers to the probability to be in state j at time t by the time we were in state i at time t − 1.

Pr{St = j | St−1 = i, St−2 = k, . . . } = Pr{St = j | St−1 = i} = pij

(4)

Now, we can define the transition probability matrix P. The latter give us all the transition probability matrix, that is to say the matrix gives us all the probability to switch from one state i to an other state j, where i, j = 1, . . . , N . P is the transition probabilities matrix. The shape of the matrix is given by the equation 5. 

p  11   p12  P= .  .  . 

p21 p22 .. .



. . . pN 1   . . . pN 2   ..  ... .  

p1N p2N . . . pN N

(5)



Note that, in our case, N = 2 because of the fact that we only have two regimes. In this particular case, we have the transition matrix given by the specification 6. 

P= 



p11 1 − p11

1 − p22   p22

(6)

Also, from this matrix, we can compute the unconditional probabilities, that to say “the probability to be in one of the two states at any given date”. We obtain the transition matrix 7.

  







π1   (1 − P22 )/(2 − P11 − P22 )   = π2 (1 − P11 )/(2 − P11 − P22 )

(7)

At this level, we would like to assess the value for all the parameters and the state.

9

2 METHODOLOGY

According to the specification of our model, we have 6 parameters to evaluate, two means and two variances parameters, plus, we have to assess the state parameters. In a first step, we will asset the state parameters. In order to do that, we use the Hamilton Filter methodology, as explain in [Hamilton(1994)]. In our example, we consider discrete state space with finite number of states. Note that, it exists three different ways to evaluate the distribution of the state.The first one is prediction, here we consider that the probability of the state depends only on the past value of rt , i.e. f (St | r1 : rt−1 ). After that, there exists a method name filtering, the latter uses all the value till time t to asset the state probability, i.e. f (St | r1 : rt ). Finally, we can use the smoothing method, here we consider all the entire training sample to determine the probability of the state variable, i.e. f (St | r1 : rT ). To assess the value of the parameter, we use the Bayesian estimation.

2.2

Bayesian estimation

Before all, it seems necessary to explain about the bayesian estimation of the parameter. Bayesians think that all parameters have a probability distribution. In fact, Bayesians condition on the data while Frequentists design procedures that work well ex ante [Nimark(2012)]. The Bayesian estimation principle is based on four principles. The first one is the likelihood principle (all the information about the parameter from an experiment is contained in the likelihood function), the second is the sufficiency principle (if a sequence of experiments is directed by a stopping rule, τ , which indicates when the experiment should stop, inference about θ should depend on τ only through the resulting sample), the third one is the conditionality principle and the fourth one is the stopping rule principle. The main components in the Bayesian inference is first the observable data, in our case the return of the electricity price or the IPI return. second, a model is required, the latter imply that we need parameters θ ∈ R+ , a prior distribution

10

2 METHODOLOGY

p(θ) : Rk → R+ . Plus, we need to write a likelihood function p(Z | θ) : RT ×n Rk →

R+ . To end, we will simulate the distribution of our parameters as explain in the following section. Of course, the first step of the simulation process is to choose a prior for the distribution of ours parameters, in other words associate to our parameters the most suitable distribution function (regarding to the range of possible value the parameter can take...). Afterwards, the Bayesian use the Bayes’ theorem based on the following equality 8:

P (θ | Z)P (Z) = P (Z | θ)P (θ) ⇐⇒ P (θ | Z) =

P (Z | θ)P (θ) ∝ P (Z | θ)P (θ) P (Z) (8)

Since P (Z) is constant (conditional on a particular model), we can use P (Z | θ)P (θ) as the posterior likelihood (a likelihood function is any function that is proportional to the probability). As explain in the next section, we will provide each parameter with a suitable distribution function.

2.3

Simulation and estimation processes

First of all, we generate a simulated dataset with known parameters. Both the probabilities, the means and the variance are known when we generate the process. The code given in the appendix 6.1.1 provide an example of the procedure. 2.3.1

Gibbs sampling methodology

As explain in the lecture notes [Lam(2008)], the Gibbs sampling method is a Markov Chain-Monte Carlo (MCMC) method. This method is used to approximate the join distribution of the parameters when direct sampling is difficult. The

11

2 METHODOLOGY

Gibbs algorithm generates a Markov chain. From the latter Markov-Chain, we would like to sample from the joint distribution. In fact, thanks to the Gibbs sampling methods, we would like to the full conditional parameter distribution for all parameters. From prior about the distribution of the parameters, we try to find the full conditional distribution of all the parameters of our model. We have to iterate the 4 steps algorithms to reach the solution. The steps are the following: 1. Write out the full posterior ignoring constant of proportionality; 2. Pick up a block of parameters θi and drop everything that does not depend on θi ; 3. Try to find the distribution Pr(θi | θ−i , y); 4. Repeat steps 2 and 3 for all parameters. Then, we obtain the full conditional distribution of all our parameters. Afterwards, we enter a loop, the Gibbs sampler steps: 1. We fix a set of starting values for all our parameters, θ(0) ; 2. Start with one of the parameters the iteration and draw a new value for the (0)

parameter from the full conditional Pr(θi | θ−i , y), θ(1) ; 3. Draw a value for an other parameter, from the actualized full conditional (0)

(1)

Pr(θj | θ−j , θi , y); 4. Iterate the previous step for all the parameters we are interesting in; 5. Repeat the steps until we get M draws, with each draw being a vector θ(t) . At the end of the day, we obtain a Markov chain of value for the parameter θ that are approximatively from our posterior.

12

2 METHODOLOGY

2.3.2

Metropolis-Hasting algorithm

In the present thesis, the Metropolis-Hasting (MH) methodology, as explain in [Robert et Casella(2004)] is used to determine the correlation coefficient between two variables. Especially, we use the methods developed by Nicholas Metropolis et al. in the copula part of our work to determine the value of the correlation coefficient between the two variables of interest. We can explain the MH as follow. In the first step, we determine the target distribution π and {X (t) }Tt=1 a sequence of random variable generated with a dependent MH with proposal q(. | X (t) . After the t-th iteration, given X (t) 1. Generate Y (t) ∼ q(y | X (t) )

2. Take X (t+1) =

  

Y (t) with probability α(X (t) , Y (t) )

 

X (t) with probability 1 − α(X (t) , Y (t) )

where α(x, y) = min

n

o

π(y)q(x|y) ,1 π(x)q(y|x)

Apply in our situation, we are in the described situation. After estimation of all the parameter inside the MS chain, we would like to assess the value of the correlation coefficient. Naturally, the latter can take value between −1 and 1. We suppose that the value for the correlation is uniformly distributed over the range [−1, 1]. At every iteration of the algorithm, a new random value generated from the uniform for the parameter will be test. Afterwards, the new value will be test against the previous tested value, inside the probability density function. We compute the value of the ratio between the likelihood of our function of interest with the new parameter and the value of the likelihood function apply with the previous value

13

2 METHODOLOGY

of the parameter. If the ratio is lower than 1, the value of the ratio become the acceptance ratio. After that, the acceptance ratio is compared to a random value obtain from the uniform 0 to 1. If the acceptance ratio is greater than the generate value, we save the generated value of our parameter of interest. If not, we save the past generated value for our parameter as the present value. We iterate that procedure, in our case the same number of time than the Gibbs algorithm. Note that, the density use for the ratio in the MH is the density of the copula model, between the two asset if interest with a bivariate normal as prior of the copula. The density of the copula is obtain by considering the first derivative of the copula according to all the parameters. Note that the Gibbs sampling and the Metropolis-Hasting algorithms are provide in annexe 6.1.2. 2.3.3

Prior distributions of the parameters

First of all, we have to choose the distribution of all our parameters, the prior distribution, regarding to the value they can take. We define the distribution of our parameters in the equation 9.

               

µ1 ∼ N (0, 100)

              

σ2 ∼ IG( 21 , 12 )

µ2 ∼ N (0, 100) σ1 ∼ IG( 21 , 12 )

(9)

P11 ∼ B( 21 , 21 ) P22 ∼ B( 12 , 21 )

where N (a, b) denotes a normal distribution with mean a and variance b, IG(a, b) and inverted gamma distribution with parameters a and b, and B(a, b) a beta distribution with parameters a and b. The equation 10 provide us with the expression of the mean parameters distribution. This equation provides us with the posterior

14

2 METHODOLOGY

full conditional distribution. Note that, for this equation, there is a identification constraint. The latter guarantee the fact that the mean of the return in the regime 1 is strictly lower than the one in the regime 2. . 

−1 1  f (µ1 , µ2 | r1 , . . . , rT , S1 , . . . , ST ) ∝ ΠTt=1 √  exp  2 2πσ(St )  1  −1 µ1 √ exp 2 10  2π10 

2 !

1  −1 µ2 √ exp 2 10  2π10 



−1 ∝ ΠTt∈τ1 exp  2 

−1 ΠTt∈τ2 exp  2

(rt − µ2 ) σ2 

(rt − µ(St )) σ(St )

(rt − µ1 ) σ1

2 !

!2 

1(µ1 < µ2 )

!2 

−1 µ1  exp 2 10

!2 



−1 µ2  exp 2 10 

2 !





2 !

1(µ1 < µ2 ) 



X −1 X 2 µ2 ∝ exp  2  rt − 2µ1 rt + T1 µ21  − 1  2σ1 t∈τ1 200 t∈τ1 







!



X −1 X 2 µ2 exp  2  rt − 2µ2 rt + T2 µ22  − 2  1(µ1 < µ2 ) 2σ2 t∈τ2 200 t∈τ2 





X rt −1  2 T1 1  ∝ exp  µ1 + − 2µ1  2 2 2 σ1 100 t∈τ1 σ1 



!





X rt −1  2 T2 1  1(µ1 < µ2 ) exp  µ2 + − 2µ2  2 2 2 σ2 100 t∈τ2 σ1 P

rt ∝ N  t∈τ21 σ1 P

rt N  t∈τ22 σ2

!

T1 1 + 2 σ1 100

!−1

T1 1 , + 2 σ1 100

T2 1 T2 1 + −1, + 2 2 σ2 100 σ2 100

!−1  

!−1 

 1(µ1

< µ2 ) (10)

Where τj = {t = 1, . . . , T | st = j}. At the end of the day, both mean parameters are normally distributed. The expression for the variance, according to the distribution of all the parameters, especially according to the actualized mean values, the density of the variance is derive from the equation 11

15

2 METHODOLOGY





rt − µj σ(st )

1 1  f (σ1 , σ2 | r1 , . . . , rT , S1 , . . . , ST ) ∝ ΠTt=1 √  exp − 2 2πσ(St )  1 σ(st )2 τ



!α+1

!2  

−β exp σ(st )2

!

! j +α+1) 1 (2 1X 1  − β + (rt − µj )2  exp 2 2 σj σj 2 t∈τj 









∝ IG α +

1 τj ,β + (rt − µj )2  2 2 t∈τj X

(11) Afterwards, it remains the probability variables, specially we would like to find the distribution of the probability according to the parameters, especially according to the actualized mean and variance values. For the probability, the prior distribution is the beta distribution.Remember that, the Beta distribution has the following expression, for α > and β > 0 we have f (x | α, β) =

1 xα−1 (1 B(α,β)

− x)β−1 1[0,1] (x),

where B is the beta distribution use to normalize the density and make integral of the latter over the period equals to 1. The equation 12 give use the form of the conditional distribution of the probability.

ξ ξ1t−1

f (p11 | r1 , . . . , rT , S1 , . . . , ST ) ∝ ΠTt=1 P111t

ξ ξ + 1 −1 t 1t 1t−1 2

P

1

−1

1

2 (1 − p11 ) 2 −1 (1 − P11 )ξ2t ξ1t−1 p11

P

1

(1 − P11 ) t ξ2t ξ1t−1 + 2 −1 X 1 1 X ξ2t ξ1t−1 + ) ∝ B( ξ1t ξ1t−1 + , 2 t 2 t

∝ P11

(12) After that, it remains the forecast of the state. We know that we would like to model discrete state space with finite number of states, in our case 2: state 1 or state 2.

16 2.3.4

2 METHODOLOGY Hamilton filter

The Hamilton filter provides us with the forecast state of the series according to the probability and the parameters. It exists several methods to draw the distribution of the state: the prediction f (St | Y1:t−1 ), the filtering f (St | Y1:t ) and the smoothing f (St | Y1:T ). We will use the latter in our study to determine the state distribution, in the smoothing technique, all the data are used to determine the distribution of the state. To determine the state value, we use the forward filtering backward sampling algorithm. As we can see in the appendix 6.1.3 and 6.1.4 is the first step of the Hamilton filter procedure. The ffbs procedure compute the log likelihood over the data and then “feed” the Hamilton algorithm. The Hamilton algorithm clearly compute the state value according to the likelihood value it receive. Note that, there is a normality assumption, to say the observation of our data are normally distributed according to the value of the parameters. Now, we have simulate and estimate all the parameters of our model, the two means, the two variances and two state probabilities (from the two probabilities, we are able to compute the two alternatives). We obtain the posterior inference for the parameters. The mean, the variance and the state parameters of our model are equal to the mean of the estimated value of the parameters estimate inside the Gibbs procedure. Note that, as you can see in the code 6.1.5, in the appendix, for the state value, we add 1 to the estimated value because of the fact that we have two states: 1 and 2.

2.4

Copula

2.4.1

Theory

Given two variables Y1 and Y2 with marginal distribution functions respectively FY1 (y1 ) and FY2 (y2 ), the copula function is way to constructing the joint distribu-

17

2 METHODOLOGY

tion of (Y1 , Y2 ). In the paper [Sklar(1959)], Sklar shown that there always exists a bivariate function C : [0, 1]2 → [0, 1] such that F (y1 , y2 ) = C(FY1 (y1 ), FY2 (y2 )), where C is also a distribution function with uniform margins on [0, 1], this function is named the copula function. The latter provide us with a bind distribution F . Note that, as mentioned in the article [Smith(2011)], ”the objective of copula modeling is not to find the function(s) C that satisfy the Sklar’s representation [...], the objective is to construct a joint distribution F from a copula function C and marginal models fo F1 and F2 “. In fact, C accounts for dependence between Y1 and Y2 . Copula model is useful because it allow us for modeling the margins separately from the dependence structure. There exist two approach is copula modeling. The first one named ”Bottom-up“. The latter consist in selecting all the univariate margin and then after dependence is introduced by an appropriate copula function C. The second approach named ”top-down“, in this situationfirst we select the joint distribution function F and after, the latter will determine the marginals. It is clear in our case that we are in the ”bottom-up“ configuration. Nelsen in the book [Nelsen(2006)] provide us with a very theoritical definition of what should be a copula function. The author provides us with a two-dimensional domain example of copulas function. Definition 2.1. A two-dimensional subcopula is a function C 0 with the following properties: 1. Dom C 0 = S1 × S2 where S1 and S2 are subsets of I containing 0 and 1; 2. C 0 is grounded and 2-increasing; 3. For every u in S1 and every v in S2 , C 0 (u, 1) = u and C 0 (1, v) = v.

18

2 METHODOLOGY

Note that for every (u, v) in DomC 0 , 0 ≤ C 0 (u, v) ≤ 1, so that Ran C 0 is a subset of I. Definition 2.2. A two-dimensional copula is a 2-subcopula C whose domain is I2 . Equivalently, a copula is a function C from I2 to Iwith the following properties: 1. For every u and v in I, C(u, 0) = 0 = C(0, v) and C(u, 1) = u and C(1, v) = v; 2. For every u1 , u2 , v1 and v2 in I such that u1 ≤ u2 and v1 ≤ v2 , C(u1 , u2 ) − C(u2 , v1 ) − C(u1 , v2 ) + C(u2 , v2 ) ≥ 0. 2.4.2

The Gaussian copula model

In the article [Pitt et al.(2006)Pitt, Chan, et Kohn], the authors suggest a method to estimate the Gaussian copula model. They write the Gaussian copula C as express in equation 13. We present the special bivariate case.

C(u | Σ) = Φ2 {Φ−1 (u1 ), Φ−1 (u2 ) | Σ}.

(13)

Where Φ is the standard normal cumulative distribution function and Φ2 (x | Σ) is the cumulative distribution function for a multivariate normal vector x having zero mean and covariance matrix Σ. Hence, the density of the Gaussian copula is given by the equation 14. 1 1 1 δ2C =| Σ |− 2 exp − (x0 Σx) exp x0 x δu1 δu2 2 2









(14)

19

2 METHODOLOGY 



 0  Where uj = Φ(xj ), x ∼ N2 (0, Σ), 0 the vector of mean, in our case   and 0 



1 ρ  Σ =    is the correlation matrix between the x variables. Indeed, uj = ρ 1 

Φ(xj ) ⇐⇒ X =  

−1



Φ (u1 )   Φ−1 (u2 )

After that, as we mentioned above, the authors use the ”bottom-up“ method. They use the marginal to obtain the the multivariate distribution. Let’s take a look at the marginal distributions. By taking yj = F −1 (uj ), for j = 1, . . . , p, where u ∼ C, they obtain the Gaussian copula regression model. In our case, because of the fact that we have a MS model, we have to consider the following expression, given by the equation 15.

yjt = F −1 (ujt ) =⇒ ujt = Fjt (yjt ) , for t = 1, . . . , T xjt = Φ−1 (ujt ) = Φ−1 (Fjt (yjt )) ⇐⇒ yjt = Fjt−1 (Φ(xjt ))

(15)

= h−1 jt (xjt ) Note that, they suppose that we have Fjt (yjt ) = Fj (yjt | sjt , θj ). The latter means that ”the marginal distribution of the j th component of case i (i = 1, . . . , n) is the same for all cases, and depends on the vector parameter θj and the m × 1 vector of covariates zij .

f (r1 , . . . , rT | µ1 , µ2 , σ1 , σ2 ) =

X s1 ∈1,2

···

X

(ΠTt=1 f (rt | µ(St ), r(St )))

(16)

sT ∈1,2

f (r1 , . . . , rT , S1 , . . . , ST | µ1 , µ2 , σ1 , σ2 ) = ΠTt=1 f (rt | µ(St ), r(St ))

(17)

20

3 DATA

Given rj,1:T = (rj,1 , . . . , rj,T ), we can write the likelihood. In our case:

T L(r1,1:T , r2,1:T , s1:T | θ) = πt=1 c(F1t (r1t ), F2t (r2t ) | Σ)f1t (r1t )f2t (r2t ) ξ

ξ

Πj=1,2 pj1jt−1 it (1 − pi1 )ξjt−1 ξit

(18)

Where ξit = 1{j} (st ). Note that ξ2t = 1 − ξ1t . Thank to the proportionality properties, the equation 19 provide use with the probability of the covariance parameter.

T c(F1t (r1t ), F2t (r2t ) | Σ) p(ρ | r11:T , r21:T , s1:T , θ−ρ ) ∝ πt=1

(19)

Note that, in our case, c will be the density for the bivariate normal distribution. The latter, as mentioned before, will be used in the MH algorithm to determined the value of the parameter ρ.

3

Data

In this section, we would like to present the data. As mentioned above, two different types of data, the Industrial Production Index (IPI) and the data of the demand of electricity in the Australia.

3.1

Industrial Production Index

The IPI is ”an economic indicator monthly released. The indicator measures the amount of output from the manufacturing, mining, electric and gas industries.“1 We obtain the IPI index for the Euro Area on the Eurostat web site2 . The data are monthly released and span the period from January 1991 to March 2014. We have 279 observations. The reference data is May 2010. 1 2

http://www.investopedia.com/ European Central bank (ECB)

21

3 DATA

For the American data, we obtain the data from the Federal Bank of Saint-Louis website3 . The reference period is 2007.

Figure 1: Industrial Production Index (IPI).

As we can observe on the graph 1, on the whole period, the trend of the two series are positively oriented. After the economic crisis, the IPI has bounced back. Note that, more than five years after the beginning of the economic crisis of 2008, both of the IPI are below their respective level of the pre-crisis period, from the beginning until May 2008. It is interesting to note that, on the two subperiods, the USA IPI was more well oriented. Indeed, if we observed the growth of the index on the two subperiods, from the beginning till the period just before the crisis the growth rates were respectively 29.7% and 62.8% On the recovery period, from June 2009 till the end, the growth rates were 8.8% and 23%. The differences between the IPI in the Euro Area and the USA are presented in the table 1. Be careful, the index are different, because of that, it does not make any sense to compare 3

Federal Reserve Bank of Saint Louis (FRED)

22

3 DATA

the mean, the minimum and the maximum between the two series. Euro Area

USA

Mean

96.6

85.9

Max.

114.7

103.3

Min.

79.5

60.2

Std.

8.5

12.0

Skew.

-0.12

-0.76

Kurt.

2.24

2.32

ADF

0.87

0.99

JB

0.03

0.00

Table 1: Descriptive statistics IPI.

As mentioned above, we would like to study the dependence of the IPI between the USA and the EA. As we know, the IPI in every country is monthly released. The graph 2 provides us with the growth rate of the IPI month on month, quarter on quarter and at the end year and year. Clearly, the volatility of the IPI growth rate is positively correlated with growth rate span. In other words, the volatility of the Month On Month (MOM) growth rate is lower than the volatility of the Quarter On Quarter (QOQ) growth rate which is lower than the volatility of the Year On Year (YOY) growth rate.

23

3 DATA

Figure 2: IPI monthly, quarterly and annually growth rate.

Finally, lets take a look at the table of the descriptive statistics of these three series. As we can observe in the table 2, the larger the time span of the considering growth rate, the more the volatility. As we can see, the distributions are stationary, we reject the null hypothesis in the Augmented Dickey-Fulller (ADF) test, plus the series are not normal, we reject the null of the Jarque-Bera (JB).

24

3 DATA

MOM

QOQ

YOY

Euro Area

USA

Euro Area

USA

Euro Area

USA

Mean

0.05

0.19

0.00

-0.01

3.03

0.08

Max.

2.24

2.06

12.60

9.60

2046.18

1776.56

Min.

-4.38

-4.30

-13.20

-15.18

-2005.48

-2007.67

Std.

0.97

0.65

4.37

2.64

776.91

504.83

Skew.

-0.86

-1.85

0.00

-0.48

-0.01

-0.01

Kurt.

5.67

12.57

2.95

7.75

2.58

4.41

ADF

0.00

0.001

0.00

0.00

0.00

0.00

JB

0.00

0.00

0.00

0.00

0.00

0.00

Table 2: Descriptive statistics of the growth rate of the IPI.

3.2

The electricity demand in Australia

We obtain the historical data of the demand of electricity on the Australian Energy Market Operator (AEMO) website 4 . The data span the period from January 2011 to April 2014. The data are released every 30 minutes. We decide to use a daily frequency of the data, we only use the observations at 7:30 am every day be cause of the fact that it is at that time that the demand is the highest, we have 1216 observations. Because of the fact that we would like to make an bivariate analysis, we decide to run the copula analysis over the data for the New South Wales and the Victoria states. Note that the data are express in megawatts (MW). The graph 3 illustrates the electricity demand in the New South Wales and in the Victoria state. As we can observed, there are saisonnalities and fluctuations of the demand of electricity. Plus, there is a negative trend of the demand in the New South Wales, at the same time the demand for electricity is stable in the 4

Australian Energy Market Operator (AEMO)

25

3 DATA

Queensland.

Figure 3: Daily demand of electricity in New South Wales and in Victoria state.

Lets take a look at the characteristics of the data sample. As we can observe in the table 3. The features of these two series are closed in term of skewness, close to 0, and kurtosis coefficient, close to 2.2. The coefficient of skewness and kurtosis must be equal, respectively to 0 and 3 by the time the data are normally distributed. Note that, both of the series are stationary, the can reject the null hypothesis at 10%.

26

3 DATA

NSW

Victoria

Mean

8134.7

5572.8

Max.

10862.1

6654.9

Min.

5663.9

4430.5

Std.

1008.9

506.7

Skew.

-0.21

-0.45

Kurt.

2.255

2.13

ADF

0.049

0.10

JB

0.00

0.00

Table 3: Descriptive statistics of the demand of electricity series.

Now, lets take a look at the daily growth rate of the electricity demand. We consider the first difference logarithmic of the daily observation, i.e. the growth rate. As we can see on the graph 4, there is not any period of high and/or high volatility. Note that we multiply the logarithmic difference by 100 to obtain directly the daily percentage change.

27

3 DATA

Figure 4: Daily growth rate of the demand of electricity in New South Wales and in Victoria state.

If we take a look at some statistics present in the table 4, we can observe that the mean of the return over the period is close to zero, which mean that the demand doesn’t change from one period to the following. Also, it is clear that the distribution of the daily growth rate is close to the normal despite the fact the Jarque-Bera test reject the null that the value are normally distributed.

28

3 DATA

NSW

Victoria

Mean

0.02

0.02

Max.

32.25

28.58

Min.

-29.37

-24.60

Std

11.94

9.89

Skew.

0.68

0.58

Kurt.

3.45

3.46

ADF

0.00

0.00

JB

0.00

0.00

Table 4: Descriptive statistics of the growth rate of the demand of electricity series.

Now, if we observe the distribution of the returns of the demand of electricity, we obtain the results presented on the graph 5. It is really more clear that the value are not normally distributed. The distribution of the return seems closer to a multivariate distribution. It seems that it is necessary to take into account 3 different distributions to fully characterized the distribution of the returns. In both of the series, the majority of the observations are concentrated around 0,

29

4 RESULTS

Figure 5: Histogram of the daily growth rate of the demand of electricity in New South Wales and in Victoria state.

4 4.1

Results Simulation process

Before all, we would like to provide the results of the simulation process. The aim of the simulation process was to give us an idea of the performance both o for the simulation process and for the parameters. By simulation process, we mean that we would like to know if the data generation process (DGP) is similar to the DGP observed in the reality, and if the algorithm succeed to find the true value we are looking for. On one hand, the graph 6 illustrate the two simulate data sets. Note that we use a 1000 iterations process for the Gibbs sampling.

30

4 RESULTS

Figure 6: Simulated data and simulated and forecast state.

As we can see on the graph 6, we arbitrary choose to define the first asset as our benchmark. The latter is distributed for the first state around −1 and for the second state around 0.1. For the variance, in both state equals to 0.1. As we can observe, the second asset is riskier than the benchmark one. The variance is bigger and the simulated value are more volatile. For the probabilities, the estimations are good in both case. The model provide us with a suitable tool to decide about the state of the economy. Plus, as we can observe, the difference between the 97.5 quantile and the 2.5 quantile is larger after the Gibbs simulation process when the variance is greater. On an other hand, we observe on the graph 7, the estimate state during the Gibbs simulation. The two states are clearly separated even if the distinction is less stable by the time the variance in the state increase.

31

4 RESULTS

Figure 7: State value beyond the Gibbs sampling.

After that, we were interesting in the convergence of the probabilities parameters. The graphs 8 provide us with the results of the estimation of the value of the parameter thanks to the Gibbs sampling simulation process.

32

4 RESULTS

Figure 8: State value beyond the Gibbs sampling.

Note that, it also possible to give the vale of the invariant probabilities. For the asset 1 we obtain π1 = 0.43 and π2 = 0.57. For the asset 2, we obtain the following results, π1 = 0.39 and π2 = 0.61. The invariants probabilities are closed. The latter means that, in both situation, the probability to be in state 1 is lower than to be in state 2. It is clear that both of the estimation of the parameter are convergent. Indeed, iteration after iteration, the cumulative mean of the estimate value of the parameter tend to be closer and closer to the truth, arbitrary chosen, value. This observation is similar with the one we can do for the mean and the variance value over the simulated samples, as illustrate thanks to the graphs 9 and 10.

33

4 RESULTS

Figure 9: Conditional mean value beyond the Gibbs sampling.

On the graph 9, it is clear that there are more uncertainty about the value of the parameter for the second asset than with the benchmark. As we can see, for the first asset, the algorithm converge very fast toward the true value. For the second asset, because of the fact that the volatility of the observations is higher, the convergence toward the true value take more time, step after step, the process meander around a mean, which is in fact the chosen value, defined at the beginning of the DGP.

34

4 RESULTS

Figure 10: Conditional variance value beyond the Gibbs sampling.

Similarly, for the variance parameter estimation process, as we can observe on the graph 10, the algorithm need more time to find the value of the parameters. Note thhat even if after 1000 iterations, in both cases the computed value meander around the mean, the true value, the fluctuation around the mean are smaller for the safer asset than for the riskier asset. Finally, we would like to compute the dependence between the two variables. At the beginning of the DGP, we define the correlation coefficient equal to 0.95. The result obtained after estimations of the parameter a thanks to the MH algorithm is observable on the graph 11.

35

4 RESULTS

Figure 11: Estimation of the correlation coefficient, copula part

As we can see, the value of the correlation coefficient converge to a value near 0.80. Note that, the performance of the model after 1000 iterations are quite good. It is clear that the spikes tend to reduce the value of the estimated parameter. The table 5 provides us with the general results from the estimations process. As we can see, the results are closed to the true value.

Asset 1

Asset 2

Probabilities

Mean

Variance

P11 =0.86 P21 = 0.11

-1.00

0.25

P12 =0.14

P22 =0.89

0.10

0.24

P11 =0.85 P21 = 0.09

-0.94

0.58

P12 =0.15 P22 = 0.91

0.06

0.44

Table 5: Convergence results.

36

4.2

4 RESULTS

On the Industrial Production Index (IPI)

If we apply our estimation procedure over observed data, especially on the MOM growth rate. A priori, because of the fact that the IPI MOM measured does not have a high volatility, one of the two states must suffer from a lack a accuracy, only few data are in one of the two regimes. Indeed, as illustrate in the graph 12, it is clear that one state dominate the other one, in our case the state 2.

Figure 12: Growth rate of the IPI in the Euro Area and in the USA, forecast state and probabilities.

As we can observed, because of the fact that the observations are more volatile in the USA, the quantile of the distribution are clearly distinct. Plus, as we can see, during the period of observation, the probability plunged but because of the fact that the probability remain greater than 0.05 there is no switch in regime. For the state of process, the graph 13 illustrates the fact that the state on the EA is stable. For the USA, the results are less clear because of the higher volatility of the data.

37

4 RESULTS

The overall results of this observations point out that the model is able to detect the crisis period, i.e the switch in the regime.

Figure 13: Estimate state of the IPI series.

Now, let’s take a look at the conditional mean and variance estimation. As mentioned above, the estimation of the parameters in the state in which the model doesn’t explore show us high volatility. Indeed, as we observed on the graphs 14 and 15, the value of the parameter converge toward a value but with a high volatility.

38

4 RESULTS

Figure 14: Estimation of the probabilities in each state for the two EA and the USA.

Figure 15: Estimation of the mean in each state for the two EA and the USA.

39

4 RESULTS

Figure 16: Estimation of the variance in each state for the two EA and the USA.

Finally, it is interesting to observed the value of the parameters at the end of the process, especially the mean value of the parameter over all the period, considering the state of the parameters. The table 6 presents the results of the convergence. It is interesting to not that, as we expected the first state, which has the lowest mean is more volatile. Probabilities

Mean

Variance

P11 =0.59 P21 = 0.02

-2.88

1.60

P12 =0.41 P22 = 0.98

0.12

1.20

USA P11 =0.66 P21 = 0.02

-1.46

1.92

P12 =0.44 P22 = 0.96

0.98

0.73

EA

Table 6: Convergence results.

In the Euro Area, the state one is a recession state with a mean about −2.88

40

4 RESULTS

and a variance around 1.60. In the second state, the mean equals to 0.12 and the variance equals to 1.20, it is a moderate expansion period. Now, lets take a look at the copula part. As we can observe on the graph 17. As observed previously on the simulated sample, the value of the correlation coefficient meander around the mean, especially the value of the correlation coefficient tend to 0.43 after 1000 iterations.

Figure 17: Estimation of the correlation coefficient, copula part.

4.3

On the Australian electricity market

As we did for the IPI study, we follow exactly the same procedure for the electricity demand analysis and will go through the output. It is clear that there are not any special situation, there is not high volatility periods followed by low volatility periods. On the graph 18, as we can see, there is no change of the state for the first state. Note that, for this estimation, we change the value of the variance parameter in

41

4 RESULTS

the mean specification. In this situation, we put the standard deviation equal to 100 instead of 10 as previously, because of the fact that the range of the fluctuation is larger. As we can see on the graph 18, it clear that the specification is well define for the two states. Plus, during all the period, there are changes of regime from 1 o 2 and 2 to 1.

Figure 18: Growth rate of the electricity demand, forecast state and probabilities.

Lets take a look at the main other results for the probabilities, the means and the variances estimation. According to the results of the state estimation, it ought be normal to observe a convergence of the parameter in the two states. For the probabilities, the graph 19 illustrates the convergence of the parameters.

42

4 RESULTS

Figure 19: Estimation of the probabilities in each state for the NSW and the Victoria states.

As we can observe, for the observations on the state of NSW, the probability to stay in the state 1 day of day is very high, around 0.9, the probability to stay in state 2 is lower, close to 0.1. In the state of Victoria, the observations about the probabilities are closed, the probabilities to stay in state 1 are high, around 0.9 and the probabilities to stay in state 2 are lower close to 0.1. For the estimation of the mean parameter, the graph 20 provide us with the results of the estimation. As we can observe,the mean parameters converges toward a value in both state. In the NSW, the mean converge toward -3 in the state 1 and around 20 in the state 2. In the state of Victoria, similarly, the means converge to a negative value in the state 1, around -3, and toward a positive value, 30, in the state 2. For the state 1, the value of the parameter converge directly after 300 iterations, in the state 2 there is convergence arrive after 300 iterations.

43

4 RESULTS

Figure 20: Estimation of the mean in each state for the two EA and the USA.

On the graph 21 illustrates the convergence of the variance parameters value. For the two series, the value of the variance converge toward a value. In the two Australian states and in both of the estimate state the variance are around 10.

44

4 RESULTS

Figure 21: Estimation of the variance in each state for the two EA and the USA.

To summaries, the table 7 provide use with the average of all the parameters over the Gibbs algorithm. As we can read from the table, in the Euro area, in the first state the mean equals to −2.15 and the variance equals to 13.9. In compare, in the second state, the means equals 16.9 and the variance equals to 9.3.In the USA, in the first state the means equals to −2.51 and the variance equals 13.2, in the second state the mean equals to 19.6 and the variance to 7.7. It is interesting to not that, as we expected the first state, which has the lowest mean is more volatile.

45

4 RESULTS

Probabilities

Mean

Variance

P11 =0.76 P21 = 0.86

-2.15

13.90

P12 =0.24 P22 = 0.14

16.97

9.35

USA P11 =0.80 P21 = 0.92

-2.51

13.18

P12 =0.20 P22 = 0.08

19.68

7.74

EA

Table 7: Convergence results.

Now lets take a look at the copula part of our study. On the graph 22, we can observe the convergence of the correlation coefficient, step after step of the MH algorithm.

Figure 22: Estimation of the correlation coefficient between, copula part.

Here, the convergence of the correlation coefficient is faster than in the simulation process. The value of the computed correlation coefficient still meander around the

46

5 CONCLUSION

mean but the volatility span is lower. There is only one spike at the beginning of the iteration process. Over the whole period, the mean of the correlation coefficient equals to 0.87. Clearly the two series are strongly positively correlated.

5

Conclusion

To conclude, as we have observed, both on the IPI study and the electricity market analysis, there is convergence of all the parameters. The latter means that we succeed to fully characterized each state of the economy in each country and each state of the electricity demand market. Plus, the copula part of the study profide us, also with a convergent estimator of the correlation coefficient. To go further, it could be interesting to assess a MS chain over the value of the correlation coefficient to observe if the value of the correlation change according to the state of the economy: the link between two economy could be different in recession period and in expansion period or by the time the demand for electricity dramatically increase or decrease. Moreover, regarding to the distribution of the demand of electricity, it could be interesting to apply directly the study over the demand instead of the growth rate of the demand. And, according to the distribution of the growth rate of the demand, it should be interesting to consider 3 states instead of 2 in the MS process.

47

REFERENCES

References [Engle et Ng(1993)] Robert F Engle et Victor K Ng : Measuring and testing the impact of news on volatility. The journal of finance, 48(5):1749–1778, 1993. [Hamilton(1994)] James Douglas Hamilton : Time series analysis, volume 2. Princeton university press Princeton, 1994. [Lam(2008)] Patrick Lam : Mcmc methods: Gibbs sampling and the metropolishasting algorithm. 2008. [Nelsen(2006)] Roger B. Nelsen : An Introduction to Copulas (Springer Series in Statistics). Springer, 2nd édition, 2006. [Nimark(2012)] Kris Nimark : Introduction to bayesian estimation. 2012. [Pitt et al.(2006)Pitt, Chan, et Kohn] Michael Pitt, David Chan et Robert Kohn : Efficient bayesian inference for gaussian copula regression models. Biometrika, 93(3):537–554, 2006. [Robert et Casella(2004)] Christian P Robert et George Casella : Monte Carlo statistical methods, volume 319. Citeseer, 2004. [Sklar(1959)] Abe Sklar : Fonctions de Répartition À N Dimensions Et Leurs Marges. Université Paris 8, 1959. [Smith(2011)] Michael S Smith : Bayesian approaches to copula modelling. 2011.

48

6 ANNEXE

6

Annexe

6.1

Matlab code

6.1.1

DGP

%%%%% Simulate DGP

mu=[-1 ;1]; sigma=[0.1;1.1]; P=[0.9 0.10; 0.10 0.90]; ip=[(1-P(2,2))/(2-P(1,1)-P(2,2));(1-P(1,1))/(2-P(1,1)-P(2,2))]; s0=1; T=200;

% number of observations

K=2;

% number of states

s=zeros(T,1); r=zeros(T,1);

i=1; s(i,1)=1+(random(’binomial’,1,P(1,2),1,1)*(s0==1)... +random(’binomial’,1,P(2,2),1,1)*(s0==2)); r(i,1)=mu(s(i,1),1)+sigma(s(i,1),1)*randn(1,1); for i=2:T s(i,1)=1+(random(’binomial’,1,P(1,2),1,1)*(s(i-1,1)==1)+... random(’binomial’,1,P(2,2),1,1)*(s(i-1,1)==2)); r(i,1)=mu(s(i,1),1)+sigma(s(i,1),1)*randn(1,1); end

%%%% set the hyperparameters

49

6 ANNEXE

a110=2; a220=2; b110=2; b220=2;

%%%% Estimation via Gibbs sampling J=10;

% number of iterations

ThetaRaw=zeros(J,6); sRaw=zeros(J,T); pfilRaw=zeros(J,T);

sRaw(1,:)=1+(r>mean(r))’; sigmaraw=0.1*ones(1,2); for j=2:J %mu muraw=zeros(1,2); muraw(1,1)=1;muraw(1,2)=0; while muraw(1,1)>muraw(1,2) T1=sum(sRaw(j-1,:)’==1); T2=sum(sRaw(j-1,:)’==2); s1=1/(T1/sigmaraw(1,1)^2+1/10); m1=s1*(sum(r.*(sRaw(j-1,:)’==1))/sigmaraw(1,1)^2); s2=1/(T2/sigmaraw(1,2)^2+1/10); m2=s2*(sum(r.*(sRaw(j-1,:)’==2))/sigmaraw(1,2)^2); muraw(1,1)=m1+s1*randn(1,1); muraw(1,2)=m2+s2*randn(1,1); end %sigma alf1=2+T1/2;

50

6 ANNEXE bet1=2+sum(((r-muraw(1,1)).^2).*(sRaw(j-1,:)’==1)); sigmaraw(1,1)=1/random(’gamma’,alf1,1/bet1,1,1); alf2=2+T2/2; bet2=2+sum(((r-muraw(1,2)).^2).*(sRaw(j-1,:)’==2)); sigmaraw(1,2)=1/random(’gamma’,alf2,1/bet2,1,1); sigmaraw=sqrt(sigmaraw);

%p11 p22 a11=a110+sum((sRaw(j-1,:)==1).*([s0 sRaw(j-1,1:(T-1))]==1)); a22=a220+sum((sRaw(j-1,:)==2).*([s0 sRaw(j-1,1:(T-1))]==2)); b11=b110+sum((sRaw(j-1,:)==2).*([s0 sRaw(j-1,1:(T-1))]==1)); b22=b220+sum((sRaw(j-1,:)==1).*([s0 sRaw(j-1,1:(T-1))]==2)); P11raw=random(’beta’,a11,b11,1,1); P22raw=random(’beta’,a22,b22,1,1); Praw=[P11raw (1-P11raw); (1-P22raw) P22raw];

ThetaRaw(j,:)=[muraw sigmaraw P11raw P22raw];

[lkl,filtp,filts,smooths] = ffbs(r,muraw’,... sigmaraw’,T,K,ip,Praw); sRaw(j,:)=smooths’; pfilRaw(j,:)=filtp(:,2)’; if round(j/10)*10==j; disp(j) end end

51

6 ANNEXE

6.1.2

Gibbs sampling

for j=2:J %%%%%%%%%%% First Asset %mu muraw1=zeros(1,2); muraw1(1,1)=1; muraw1(1,2)=0; while muraw1(1,1)>muraw1(1,2) T11=sum(sRaw1(j-1,:)’==1); T21=sum(sRaw1(j-1,:)’==2); s11=1/(T11/sigmaraw1(1,1)^2+1/10); m11=s11*(sum(r1.*(sRaw1(j-1,:)’==1))/sigmaraw1(1,1)^2); s21=1/(T21/sigmaraw1(1,2)^2+1/10); m21=s21*(sum(r1.*(sRaw1(j-1,:)’==2))/sigmaraw1(1,2)^2); muraw1(1,1)=m11+s11*randn(1,1); muraw1(1,2)=m21+s21*randn(1,1); end %sigma alf11=alf+T11/2; bet11=bet+sum(((r1-muraw1(1,1)).^2).*(sRaw1(j-1,:)’==1)); sigmaraw1(1,1)=1/random(’gamma’,alf11,1/bet11,1,1); alf21=alf+T21/2; bet21=bet+sum(((r1-muraw1(1,2)).^2).*(sRaw1(j-1,:)’==2)); sigmaraw1(1,2)=1/random(’gamma’,alf21,1/bet21,1,1); sigmaraw1=sqrt(sigmaraw1);

%p11 p22 a111=a110+sum((sRaw1(j-1,:)==1).*([s01 sRaw1(j-1,1:(T-1))]==1));

52

6 ANNEXE a221=a220+sum((sRaw1(j-1,:)==2).*([s01 sRaw1(j-1,1:(T-1))]==2)); b111=b110+sum((sRaw1(j-1,:)==2).*([s01 sRaw1(j-1,1:(T-1))]==1)); b221=b220+sum((sRaw1(j-1,:)==1).*([s01 sRaw1(j-1,1:(T-1))]==2)); P11raw1=random(’beta’,a111,b111,1,1); P22raw1=random(’beta’,a221,b221,1,1); Praw1=[P11raw1 (1-P11raw1); (1-P22raw1) P22raw1];

ThetaRaw1(j,:)=[muraw1 sigmaraw1 P11raw1 P22raw1];

ip1=[(1-Praw1(2,2))/(2-Praw1(1,1)-Praw1(2,2));(1-Praw1(1,1))/... (2-Praw1(1,1)-Praw1(2,2))];

[lkl1,filtp1,filts1,smooths1] = ffbs(r1,muraw1’,... sigmaraw1’,T,K,ip1,Praw1); sRaw1(j,:)=smooths1’; pfilRaw1(j,:)=filtp1(:,2)’;

%%%%%%%%%%%% Second Asset %mu muraw2=zeros(1,2); muraw2(1,1)=1; muraw2(1,2)=0; while muraw2(1,1)>muraw2(1,2) T12=sum(sRaw2(j-1,:)’==1); T22=sum(sRaw2(j-1,:)’==2); s12=1/(T12/sigmaraw2(1,1)^2+1/10); m12=s12*(sum(r2.*(sRaw2(j-1,:)’==1))/sigmaraw2(1,1)^2); s22=1/(T22/sigmaraw2(1,2)^2+1/10); m22=s22*(sum(r2.*(sRaw2(j-1,:)’==2))/sigmaraw2(1,2)^2);

53

6 ANNEXE muraw2(1,1)=m12+s12*randn(1,1); muraw2(1,2)=m22+s22*randn(1,1); end %sigma alf12=alf+T12/2; bet12=bet+sum(((r2-muraw2(1,1)).^2).*(sRaw2(j-1,:)’==1)); sigmaraw2(1,1)=1/random(’gamma’,alf12,1/bet12,1,1); alf22=alf+T22/2; bet22=bet+sum(((r2-muraw2(1,2)).^2).*(sRaw2(j-1,:)’==2)); sigmaraw2(1,2)=1/random(’gamma’,alf22,1/bet22,1,1); sigmaraw2=sqrt(sigmaraw2);

%p11 p22 a112=a110+sum((sRaw2(j-1,:)==1).*([s02 sRaw2(j-1,1:(T-1))]==1));... % State 1 to state 1 a222=a220+sum((sRaw2(j-1,:)==2).*([s02 sRaw2(j-1,1:(T-1))]==2));... % State 2 to state 2 b112=b110+sum((sRaw2(j-1,:)==2).*([s02 sRaw2(j-1,1:(T-1))]==1));... % State 1 to state 2 b222=b220+sum((sRaw2(j-1,:)==1).*([s02 sRaw2(j-1,1:(T-1))]==2));... % State 2 to state 1 P11raw2=random(’beta’,a112,b112,1,1); P22raw2=random(’beta’,a222,b222,1,1); Praw2=[P11raw2 (1-P11raw2); (1-P22raw2) P22raw2];

ThetaRaw2(j,:)=[muraw2 sigmaraw2 P11raw2 P22raw2];

ip2=[(1-Praw2(2,2))/(2-Praw2(1,1)-Praw2(2,2));(1-Praw2(1,1))/... (2-Praw2(1,1)-Praw2(2,2))];

54

6 ANNEXE

[lkl2,filtp2,filts2,smooths2] = ffbs(r2,muraw2’,... sigmaraw2’,T,K,ip2,Praw2); sRaw2(j,:)=smooths2’; pfilRaw2(j,:)=filtp2(:,2)’;

%%%%%%%%%%%%%% Copula part mu1t=muraw1(1,sRaw1(j,:)’)’; sigma1t=sigmaraw1(1,sRaw1(j,:)’)’; u1=cdf(’normal’,r1,mu1t,sigma1t);

mu2t=muraw2(1,sRaw2(j,:)’)’; sigma2t=sigmaraw2(1,sRaw2(j,:)’)’; u2=cdf(’normal’,r2,mu2t,sigma2t);

x=[icdf(’normal’,u1,0,1) icdf(’normal’,u2,0,1)];

%% The MH algorithm: apply at the same time than the Gibbs sampling procedure. %

We apply the same number of iteration of the process

nMH=2; myRhoMH(j,1)=myRhoMH(j-1,1); for i=1:nMH, myRhoStar=1-2*rand(1,1); mySigmaOld=[1, myRhoMH(j,1); myRhoMH(j,1), 1]; mySigmaNew=[1, myRhoStar; myRhoStar, 1]; Fnew = mvnpdf([x(:, 1) x(:, 2)],[0,0], mySigmaNew);

55

6 ANNEXE Fold = mvnpdf([x(:, 1) x(:, 2)],[0,0], mySigmaOld);

alfaxy=min([exp(sum(log(Fnew)-log(Fold))),1]);

Acc(i,1)=alfaxy; u=rand(1,1); if u0.5);

Smile Life

When life gives you a hundred reasons to cry, show life that you have a thousand reasons to smile

Get in touch

© Copyright 2015 - 2024 PDFFOX.COM - All rights reserved.