switching volatility models Bayesian changepoint and time-varying parameter learning in regime Digital Signal Processing

Contents lists available atScienceDirect

Digital Signal Processing

www.elsevier.com/locate/dsp

Bayesian changepoint and time-varying parameter learning in regime switching volatility models

M. Serdar Yümlü

^a,∗

, Fikret S. Gürgen

^a

, A. Taylan Cemgil

^a

, Nesrin Okay

^b

aDepartmentofComputerEngineering,Bo˘gaziçiUniversity,34342Bebek,Istanbul,Turkey bDepartmentofManagement,Bo˘gaziçiUniversity,34342Bebek,Istanbul,Turkey

a r t i c l e i n f o a b s t ra c t

Articlehistory:

Availableonline20February2015

Keywords:

MultipleChangepointDetection(MCD) SequentialMonteCarlo(SMC)methods ParticleFiltering(PF)

AuxiliaryParticleFiltering(APF) ExponentialGeneralizedAutoregressive ConditionalHeteroskedasticity(EGARCH) Volatilitymodeling

This paper proposes a combined state and piecewise time-varying parameter learning technique in regimeswitchingvolatilitymodelsusingmultiplechangepointdetection.ThisapproachisaSequential Monte Carlo method for estimating GARCH & EGARCH based volatility models with an unknown numberofchangepoints.Modernauxiliaryparticlefilteringtechniquesareusedtocalculatetheposterior densities andonlineforecasts.Thisapproachalsoautomaticallydealswiththecommonancestralpath dependenceproblemfacedinthesetypevolatilitymodels.ThemodelistestedonBorsaIstanbul(BIST) formerlyknown asIstanbulStock Exchange(ISE)marketdata usingdailylogreturns.Afull structural changepointspecificationisdefinedinwhichallparametersoftheconditionalvarianceofthevolatility models are dynamic. Finally, it is shown with simulation experiments that the proposed approach partitions the seriesintoseveralregimes and learnsthe parametersof eachregime’svolatility model inparallelwiththemultiplechangepointdetectionprocess.

©^{2015ElsevierInc.}All rights reserved.

1. Introduction

Thefinancialmarket isacomplex,evolving,andnonlineardy- namicalsystem.Thefinancialtimeseriesareinherentlynoisy,non- stationary,andchaotic[1].Thisbringsthealteringstructureofthe distributionoffinancialtimeseriesovertime.Themeanandvari- anceofthetimeseriesare non-stationaryinthemselves.Eventhe relationship betweendifferenttime series such ascurrency rates andstockexchangesmaychangeovertime.Financial TimeSeries Prediction (FTSP)is firstinvestigating theunderlying structureof thetimeseries,fittingamodelbasedonthedatachangingintime andusingthismodelforpredictionofthefollowingtimesteps.

Modelingsuchdynamicalandnon-stationarytimeseriesisex- pected to be a challenging task. FTSP is also a difficult signal processingproblemthathashiddenvariablesandlacksobservable datafordeterminingtheunderlyingstructureoftheseries,ifone exists.Weusedconditionalvariance(volatility)inthispaperwhich is the time-dependent heteroskedastic variance. Briefly, volatility isthe measureoftime domainchangeabilityof assetreturns [2].

Stockpricesfluctuateduetomanydifferentriskinessinthemarket causingheteroskedasticityintheobservations.Asaresult,accurate predictioncanonlybeachievedwithaforecastofthetimedepen-

*

Correspondingauthor.

E-mailaddress:serdaryumlu@gmail.com(M.S. Yümlü).

dentvarianceofthestock.Volatilitiesofassetreturnsareusedto analyzemarketrisks,portfolioselectionandmarkettiming[3].

Recent financial time series predictionstudies were based on applications ofArtificial NeuralNetworks,Mixture ofExperts[3], and Support Vector Regression [4]. In contrast, our focus is on Sequential Monte Carlo (SMC) methods to perform online pre- diction based on joint probability distribution in Hidden Markov Model (HMM) for non-linearity and non-Gaussian scenarios [5].

SMCmethods,ageneralclassofMonteCarlomethods,aremostly used for sampling from sequences of distributions. Simple ex- amples of these algorithms are extensively used in the tracking andsignalprocessingliterature.Recentdevelopmentsindicatethat these techniques have much more general applicability, and can be applied to statistical inference problems in signal processing.

Financial time seriesinclude severalswitchingregimes whichare notpossibletobemodeledbyasinglemodel.Forthisreason,we aim toaddress challengesin estimatingofvolatility modelssuch asGARCH&EGARCHmodelsthataresubjecttostructuralregime switchesbyusingaBayesianSMCmethodsapproach.

For mostreal-world problems,the optimal Bayesian inference usingrecursionsisintractableandapproximatesolutions mustbe used. Within the space of approximate solutions, the extended Kalmanfilter(EKF)hasbecomeoneofthemostwidelyusedalgo- rithmswithapplications instate, parameterand dualestimation.

Van der Merwe extended thesealgorithms to a family of Sigma Point Kalman Filters (SPKF) in dynamic state-space models [6].

http://dx.doi.org/10.1016/j.dsp.2015.02.001 1051-2004/©^{2015ElsevierInc.}All rights reserved.

Matsumoto attempt to construct a model from time-series data andmakeanonlinepredictionwhenthelinearassumptionisnot valid.TheproblemisformulatedwithinaBayesianframeworkim- plementedbytheSMCmethod[7].

Financial time series include several switching regimes, for this reasonwe have to detect the unknown number of change- points together withthe regime parameters. Previous studies on changepoint detection have the roots of Chib’s approach. Chib (1998) proposed a convenient state-space formulation of struc- turalbreakmodelsinwhichadiscrete latentstatevariabledrives structuralbreaks andfollowsa constrainedMarkov chain.An ef- ficientMarkov Chain Monte Carlo(MCMC) algorithm forestima- tion and computing Bayes factors for the purpose of model se- lection is provided in his paper [8]. Chopin (2007) has devel- oped aparticle filteringalgorithm forestimatingstructuralbreak models in which the fixed model parameters are formulated as part of the latent variables [9]. He and Maheu proposed a SMC method for estimating GARCH models subject to an unknown number of structural breaks. They applied their model to daily NASDAQ returns and followed Chib’s formulation of structural breaks[10].Turneret al.proposedanadaptivesequentialapproach forBayesianChangePointDetectionusingBayesianOnlineChange- pointDetection (BOCPD)ofMackayby introducingan Underlying PredictiveModel(UPM)andahazardfunction[11].

Tsay et al. reviewed SMC methods (particle filters (PF)) with special emphasis on its potential applications in financial time.

They analyzed Liu–West filter (2001), Storvik filter (2002) and particle learning (PL) of Carvalho, Johannes, Lopes and Polson (2010) and the Auxiliary Particle Filter (APF) of Pitt and Shep- hard (1999) [12]. Eckley and Fearnhead also analyzed different typeofchangepointmodelsintimeseriesproblemsincludingsin- glechangepoint modelsandmultiplechangepointmodels.They compared Bayesian based and other approaches in change point modeling[13].Whiteleyet al.provided explanations aboutparti- clefiltering, sequential importance sampling andauxiliary parti- cle filtering techniques. He also introduced recent developments in auxiliary particle filtering with applications to several prob- lems [14]. In empirical Bayesian change point detection, Paquet builtanonlinealgorithm ofAdamsandMackay(2007), whocast theproductpartitionintoaBayesiangraphicalmodel.Adams and Mackaycomputedtheprobability distributionofthelengthofthe current “run” or time since the last change point, using a sim- plemessage-passingalgorithm[15].Thismodelenablesaninfinite number of hiddenstates and change points in observations and withthis property it hasan advantage over Chib’ s HMM based approach where the number of change points have to specified previously[16].

Aparticle filter is a class of SMCmethods that approximates posteriordistributionofthelatentstate variablesbya setofpar- ticlesandassociated importance weights. In orderto understand the underlying switching models we propose a combined mul- tiple changepoint detection and time-varyingparameter learning in regime switching state space models. Our approach includes both state estimation and parameter estimation steps. One ap- proachto learning about fixed parameters is the mixture kernel smoothingmethod ofLiu andWest(2001). This paperaddresses thechallengesinestimatingtheunknownnumberofchangepoints involatilitymodelsusingBayesianinferenceandSMCmethodsap- proaches.Aparticlefilteringapproach basedonAPFtosequential estimationisbuiltontopofthechangepointmodelofChib(1998).

Theunknownnumberofchangepointsandtheregimeparameters areestimatedjointly.

The objective of thisresearch is to investigate the use of re- centSMCmethodsapproachesfortheMultipleChangepointDetec- tion(MCD)problemsin financialtime seriespredictionproblems without knowing the number of changepoints [17–19]. Financial

time series,esp.stock market timeseriesshow non-linearityand stochasticity. In order to solve this, several volatility models are usedincludingGARCH&EGARCHanditsderivatives.Theproposed approachisaSMCmethodthatcombinesAuxiliaryParticleFilter- ing (APF)withMixtureKernelSmoothing(MKS) fortime-varying parameterandchangepointestimationinregimeswitchingGARCH

&EGARCHbasedvolatilitymodels.RecentAPFtechniquesareused tocalculatetheposteriordensitiesandforecastsinreal-time.This approach alsoautomaticallydeals withthe commonpath depen- denceproblemofthesetypevolatilitymodels.

The contributions in this paper which shows a different so- lution to the problem of multiple changepoint detection, time- varying parameter learning and volatility modeling by extending previousapproaches[9,10,12]are:

• ^{Providing multiple regime switching state space models in-} stead of fitting a global andsingle GARCH & EGARCH based modeltothetimeseries

• ^{Combined multiple} changepoint detection and time-varying parameterlearninginregimeswitchingstatespacemodelsfor volatilitymodeling

• ^{Appliedextensive} simulationsusing GARCH& EGARCH based models,Student-tdistributionoversyntheticdatasetandvali- datedusinganemergingmarketdataofBorsaIstanbul(BIST)

With these features, the proposed approach has advantages over traditional andsingle modelapproaches [20]. The proposed solution is first testedover synthetic dataset andthen validated usingan emerging marketdataof BorsaIstanbul wherevolatility showsmorefluctuationsthanUSmarkets.BesidesGaussiandistri- butions,Student-tdistributionwhichismoresuitableforfinancial timeseriesdataisconsideredwithbothGARCHandEGARCHbased volatilitymodels.Thispaperprovidesadifferentandmulti-regime auxiliaryparticlefilteringbasedpiecewisetimevaryingparameter learningapproachinregimeswitchingvolatilitymodels.

The rest of the paper is organized as follows. In Sections 2 and3,financialtime series,volatilitymodelsandstate spacerep- resentation for multiple changepoint models are presented. SMC approachesandAPFmodelsareintroducedinSection4.Sections5 and6presentsmixturekernelsmoothingbasedparameterlearning and the proposed approach of APF based MCD modelin volatil- itymodels.Section7showsexperimentsandresultsonsynthetic data sets and real time series data of Borsa Istanbul (BIST) 100 indexand concludingremarks are givenwith severaldiscussions in the last section. The originalityof this paperlies on the pro- posedapproachthatcombinesAPFbasedSMCmethodsandKernel Smoothing based parameter learning to provide combined state andpiecewisetime-varyingparameter learninginregime switch- ing volatility models using multiple changepoint detection tech- niques.

2. Financialtimeseries&volatilitymodels

The financial time series data set is used from an emerging market,BorsaIstanbulformerlyknownasIstanbulStockExchange (ISE),BISTNational100index.BIST100reflectsallthecharacteris- ticsofthemarketbecauseitincludesthemostactiveandvolumet- ric100stockswhichare selectedamongthestocksofcompanies traded on the national market. ISE was established in 1988 and allthe datastarting from1988up to2013isgatheredtoanalyze thebehavioroftheproposedapproach.Thedataset includes6254 dailyobservationsstartingfrom04January1988to01March2013 (Fig. 1).

Turkish stock market BIST, experiencing powerful fluctuations translatedbylargepricemovements,isalightstaramongvolatile emerging markets. Closer look to the BIST reveal highdegree of

Fig. 1. BIST 100 index USD based close prices & BIST 100 log returns.

persistence and strong time dependence inconditional variances thatmakesforecastingvolatilityanimportantissue[2].

Conditionalvariances are knownto be unobservable time de- pendentfeaturesandcanbeanalyzedusingAutoregressiveCondi- tional Heteroskedasticity (ARCH) processes based heteroskedastic models proposed by Engle (1982). Bollerslev (1986) generalized this approach offering a Generalized Autoregressive Conditional Heteroskedasticity (GARCH)modelin whichconditional variances aregovernedbyalinearautoregressiveprocessofpastsquaredre- turns andvariances [2].In the following section, we will survey somecommonvolatilitymodels.

2.1. GeneralizedAutoregressiveConditionalHeteroskedasticitymodels

GARCHisamechanismthat includespastvariancesinthe ex- planation of future variances. More specifically, GARCH is a time seriesmodelingtechnique that usespastvariancesandpastvari- ance forecaststo forecast futurevariances.A GARCH(Generalized AutoRegressiveConditionalHeteroskedasticity)process isa com- mon model used in time series analysis for analyzing stochastic volatility[15].Specifically,aGARCH(1,1)modelisthefollowing:

y_k

= σ

_kz_k (1)

z_k

∼

^N

(

0

,

1

)

(2)

σ

_k²

= ω + α

y_k²₋₁

+ β σ

_k²₋₁ (3) In this model y_k is the demeaned log return series,

σ

_k² rep- resentsthevolatilityand

ω

,

α

andβ arethevolatilityparameters whichcontroltheeffectsofthelogreturnsandhistoricalvolatility.

ARCH parameter

α

and GARCH parameter β controls the persis- tenceofa GARCHmodel whichshowshowfast volatilitiesdecay aftercriticalchanges.Thismodelassumesthat

ω

,

α

andβ >0 and

α

_{+ β <}1.

GARCH hasgained a lotof interest andis widely accepted.It takesintoaccountexcesskurtosis(i.e.fattailbehavior)andvolatil- ityclustering,twoimportantcharacteristicsoffinancialtimeseries.

Itprovides accurateforecastsofvariancesandcovarianceofasset returnsthroughitsabilitytomodeltimevaryingconditionalvari- ances. As a consequence, you can apply GARCH models to such diversefieldsasrisk management,portfoliomanagement andas- set allocation,optionpricing,foreignexchange,andthetermstruc- ture of interest rates. Although GARCH models are useful across awide rangeofapplications,they dohavesome limitationssuch astheirinabilitytocaptureirregularmarketmovementsandtheir

parametric specificationthat operates better understable market conditions. GARCH models oftenfail to fullycapture thefat tails observed in asset return series.Heteroskedasticity explainssome ofthefattailbehavior,buttypicallynotallofit.

2.2. ExponentialGARCH(EGARCH)models

Asymmetric leverage volatility models incorporates the hy- pothesis that negative shocks cause morevolatility than positive shocks. News Effect causing asymmetric volatility is a necessary result to be included in the estimations [20]. Despite their suc- cessfulapplicationsARCHandGARCHmodelscannotcapturesome importantfacts inthedata.The mostimportantfact isthelever- ageorasymmetriceffectdiscoveredbyBlackandconfirmedbythe findingsofNelson[21].Thiseffectclaimsthatabadnews,adrop in priceincreasesvolatility more thanan unexpected increase in priceofsimilarmagnitude.Becausewetakesquareoftheinnova- tions inARCHandGARCH, weare unableto seizethiseffectand understand if there will be a difference from the point of view of volatility.This effectsuggeststhat a symmetryconstraint over y_k is not appropriate. Nelson proposed a model to capture such asymmetric effects called Exponential Generalized Autoregressive Conditional Heteroskedasticity (EGARCH) [21]. This model is de- finedasbelow:

log

( σ

_k²

) = ω + β

^log

( σ

_k²₋₁

) + γ

y_k−1

/ σ

k−¹

+ α (|

^yk−¹

|/ σ

k−¹

− 

2

/ π )

(4)

where

ω

,

α

,β and

γ

areconstantparameters.Thismodelisasym- metricanditisabletocovertheeffectofthesignofthereturns because yk−¹/

σ

k−¹ has a coefficient of

γ

. When this coefficient is negative, it will generate more volatility than positive return shocksbecauseofthesigneffect.TheEGARCHmodelallowsgood news and bad news to have a different impact on volatility on the other hand GARCH doesnot. Bad newshas a greater impact thangoodnews.TheEGARCHmodelalsoallowsbignewstohave a greaterimpactonvolatilitythanGARCH model.Thetime series is governed by multiple regimes andour aim is to estimate the changepoint locations and regime parameters between switching regimes.

Fig. 2showsthechangepoints,returnseriesandrealizationsby usingEGARCHmodeltoexplaineffectofbothGARCHandEGARCH parameters. The series is divided into 5 different regimes with changepointsateach100thtimestep. The1st,3rdand5thregimes

Fig. 2. Volatility realization by EGARCH parameters (changepoints, returns, volatility respectively).

have constant parameters of

ω

=^0,

α

=⁰.05 and β =⁰.05. As a result of this volatility shows smooth persistence. In the 2nd regimewefirstseetheeffectofvolatilityconstantbychanging

ω

to

ω

=¹.75.Thiscontrolsrangeofthereturns butdoesnothave anychangingeffectonthevolatility. Toseetheeffectofparame- ter

α

,theparametervalueisaset to0.85inregime 3.Parameter

α

istheARCH parameterandshowsthe effectofthelogreturns on volatility. The

α

parameter represents a magnitude effect or thesymmetriceffectofthemodel.β measures thepersistencein conditionalvolatilityirrespective ofanythinghappeningin thefi- nancialmarket.When β isrelatively large,thenvolatility takesa longtimetodieoutfollowingacrisisinthemarket.Herethehis- toricalvolatility is not incorporated yetas β is still the same as previousregimes.Inregime 4,werevertbacktheARCHparameter

α

toits original stateandincreasing historicalvolatility effectby changingtheGARCHparameterβ to0.85.Bytheeffectofhistorical volatility, thevolatility started toincrease withlittle fluctuations comingfromlogreturnswhicharecontrolledby

α

.The

γ

param- etermeasurestheasymmetryortheleverageeffect,theparameter of importance so that the EGARCH model allows for testing of asymmetries. If

γ

=^{0 then the model is symmetric (Regimes1,} 3 and 5). When

γ

<0, then positive shocks (good news) gen- erateless volatility than negative shocks (bad news) (Regime 4).

When

γ

>0,itimpliesthat positiveinnovations are moredesta- bilizing than negative innovations (Regime 2). We assume

σ

₀² at time k=^{0 is to be the} unconditional variance

ω

/(1−

α

− β)^, (y_k∼^N(0,

σ

_k²)).

3. Statespacerepresentationformultiplechangepointmodels

StateSpaceModel(SSM)ofatime seriesisusuallydefinedas s_k and y_k which are the latent state variables and the observa- tions respectively.Wewillhereusesk asthelatentstate variable which represents the regime identifier. We here define the den- sityofthe observationsas p(yk|^y1:^k−¹,θk,sk)where θk isthe set of parameters at time k. The idea behind the model is the for- mulation of the changepoint modelin terms of a latent variable thatindicates theregime fromwhichaparticularobservationhas been drawn. The regimes of the process is defined first andthe parameters θ_k are drawnby an unknown state variable s_k where sk=^m,m∈(¹,2,. . . ,M)andθk= θm whichmeanstheparameters θk ateachtimek are equalandsameasthemth regime θm.The latentvariableisadiscrete Markovprocesss_k withthetransition probabilities P limited so that the state variable can either stay inthe sameregime orjumptoa newregime [6].In speechpro- cessing,suchmodelsarecalledLeft–RightHMMs.Left–RightHMM which is also known asBakis HMM can model signalschanging overtimelikespeechsignals[22].Fig. 3showsa4-stateleft–right HMMmodel.

P

=

⎡

⎢ ⎢

⎣

p₁₁ 1

−

^p11 0 0

0 p22 1

−

^p22 0

0 0 p₃₃ 1

−

^p33

0 0 0 1

⎤

⎥ ⎥

⎦

⁽⁵⁾

This changepoint model is selected in this paper in order to give alimittotheparameterspacebutthereare alsoothertech- niquesforstatespacerepresentationofmultiplechangepointmod-

Fig. 3. A 4-state left–right HMM model.

Fig. 4. Graphical model for changepoints.

Fig. 5. Run length indicator rk[15].

els in the literature. Fearnhead assigned changepoint flags and usedthemodel representationgivenin Fig. 4.Here s_k represents the changepoint flag,

σ

k is the latent variable and yk is the ob- servation[19].[^sk=⁰]^{denotes theIversonbracketwhichactsas} an indicator functionandresultsto 1ifthe statement insidethe bracketissatisfied,and0otherwise.

s_k

∼

^p

(

s_k

|

^sk−1

),

Changepoint flags

∈ {

⁰

,

1

}

σ

_k

_{∼ [}

s_k

=

⁰

]

^f

( σ

_k

_| σ

_k−1

) + [

^sk

=

¹

] π ( σ

_k

), π ( σ

_k

)

: Reinitialization y_k

∼

^p

(

y_k

| σ

_k

),

Observations

AdamsandMackayusedcurrent“runlength”orthetimesince thelastchangepointtoestimatetheposteriordistribution.Hererk denotesthelength ofthecurrentrunattime k.Fig. 5 showsthe run lengthr_k asafunction oftime andr_k drops tozero whena changepointoccurs[11,15].

3.1. Multiplechangepointsforvolatilitymodels

Weusedpartialandfullchangepointmodelswhereonlysubset oforalltheparametersaresubjecttochange.BothGARCHandex- ponentialGARCHmodelsareusedasthevolatilityprocessesofthe timeseries.WefirstusedNormaldistributionbutempiricalstudies often suggest fat tails in the distribution of asset returns, there- fore,we extendedthe proposed solutionto incorporatestudent-t distributionwith

ν

degreesoffreedomaswell.Firstwecalculated volatility usingGARCH [15] andNelson’s EGARCH [21]. Residuals and volatility are used in the process model and measurement model equations.Measurement and process modelequations are adaptedtoGARCH andEGARCHmodels. SMCmethods areimple- mented using the measurement andprocess model equationsin

Fig. 6. Graphical representation for the multiple changepoint model.

Sections2.1and2.2.Fig. 6showsthegraphicalrepresentationfor themultiplechangepointmodeldefinedinthispaper.

Inthisgraphicalmodelskandykrepresentsthelatentvariables andtheobservationsrespectivelywheres_kistheregimeidentifier.

θ= [θ1,θ2,. . . ,θ_M¯] ^{is the set of parameters at differentregimes.}

Probabilitydensitiesofstateandobservationsaregivenbelow:

s_k

∼

^p

(

s_k

|

^sk−¹

, θ )

(6)

y_k

∼

^p

(

y_k

|

^sk

, θ )

(7)

ThestatemodelisaMarkovprocesswhereitisconditionalon the previous state and the parameters θ, and the observationat time k is conditionalon state sk.In orderto determinethe state atanygivenpoint k, itrequiresto estimatethe filteredposterior distribution p(s_k|^y1:k,θ ) recursively intime. The posterior distri- bution of the filtered states can be derived by using Bayes Rule as

p

(

s_k

, θ |

y_1:k

) =

^p

(

y_k

|

^sk

, θ )

p

(

s_k

, θ |

^y1:k−1

)

p

(

y_k

|

y_1:k−1

, θ )

⁽⁸⁾

where

p

(

sk

, θ |

y1:^k−¹

) =

p

(

sk

|

sk−¹

, θ )

p

(

sk−¹

, θ|

y1:^k−¹

)

dsk−¹ (9) and

p

(

y_k

|

^y1:^k−¹

, θ ) =

p

(

y_k

|

^sk

, θ )

p

(

s_k

, θ|

^y1:^k−¹

)

ds_k (10) WemodeltheobservationsasGARCHorEGARCHprocessesby using

ω

,

α

, β and

γ

parameters. Let θm= [

ω

m,

α

m,βm,Pmm] ^or θm= [

ω

m,

α

m,βm,

γ

m,P_mm]^{bethemodelparametersands}misthe latentvariableinstatem andθ= [θ1,θ2,. . . ,θ_M¯]^{istheparameter} spaceinall statesandregimes.sk=^sm attime stepk wherekth time step resides in regime m. The likelihood of y_k for the full changepointmodelisthendefinedas

p

(

y_k

|

^y1:^k−¹

,

s_k

, θ ) =

^p

(

y_k

|

^y1:^k−¹

, θ

s_k

)

= (

²

π σ

_k²

)

^−1/2exp

( −

^y

2 k

2

σ

_k²

⁾

⁽¹¹⁾

where yk=

σ

kzk,zk∼^N(0,1).

In full changepoint model all the parameters of the GARCH modelaresubjecttochangeasgivenbelow.

yk

= σ

kzk

, σ

_k²

=

ws_k

+ α

s_ky²_k−1

+ β

s_k

σ

_k²₋₁ (12)

θ

m

= [ ω

m

, α

m

, β

m

,

P_mm

]

⁽¹³⁾

ForEGARCHmodels:

log

( σ

_k²

) = ω

s_k

+ β

s_klog

( σ

_k²₋₁

) + γ

s_kyk−¹

/ σ

k−¹

+ α

s_k

(|

yk−¹

|/ σ

k−¹

− 

2

/ π )

(14)

andθm= [

ω

m,

α

m,βm,

γ

m,Pmm]^.

Inpartialchangepointmodelonly parameter wsk issubjectto change.GARCHmodelisthendefinedasfollows:

σ

_k²

= ω

s_k

+ α

y²_k−1

+ β σ

_k²₋₁ (15)

θ

m

= [ ω

m

, α , β,

Pmm

]

(16)

whereyk=

σ

kzk,

σ

_k²=

ω

s_k+

α

y²_k−1+ β

σ

_k²₋₁. ForEGARCHmodels:

log

( σ

_k²

) = ω

s_k

+ β

^log

( σ

_k²₋₁

) + γ

y_k−1

/ σ

_k−1

+ α ( |

^yk−1

|/ σ

_k−1

₋ ^

2

/ π )

(17)

andθm= [

ω

m,

α

,β,

γ

,P_mm]^.

Aswementionedaboveempiricalstudiesoftensuggestfattails inthedistributionofassetreturns,therefore,weextendedthepro- posedsolutiontoincorporatestudent-tdistributionwith

ν

degrees offreedom.

Letθm= [

ω

m,

α

,β,

ν

,Pmm]orθm= [

ω

m,

α

,β,

ν

,

γ

,Pmm]bethe modelparametersandsmisthelatentvariableinstatem andθ= [θ1,θ2,. . . ,θ_M¯] ^{is the parameter space inall states and regimes.}

Thedensity of y_k for thepartial changepoint modelwhere y_k=

σ

kz_k,z_k∼^N(0,1)is

p

(

yk

|

y1:^k−¹

,

sk

, θ ) =

p

(

yk

|

y1:^k−¹

, θ

s_k

)

= (( ν +

¹

)/

2

) σ

k

(θ/

2

) √

π ν ⁽

¹

⁺

y_k²

νσ

_k²

⁾

−

( ν +

¹

2

)

(18) wherey_k=

σ

kz_k,

σ

_k²=

ω

s_k+

α

y²_k−1+ β

σ

_k²₋₁.

Full changepoint model and EGARCH models applies as the sameinthenormaldistribution.

4. SequentialMonteCarlo(SMC)methods

Most of the real world problems require non-linear, non- Gaussian scenarios and in most real-world applications, the in- tegrations with respect to sk−¹ and sk in (9), (10) and the im- plementation of Bayes’ theorem in (8) are both analytically in- tractable and/or computationally costly. Hidden Markov Models (HMM)andLinearDynamicalSystems(LDS)systemsfailintopro- vidingtractablefilteringsolutionsforthesetypeofnon-linearand non-Gaussian scenarios. In this paper, a probabilistic time series modeling and changepoint detection approach for stock market time seriesusing non-linearand non-Gaussian noise scenarios is introduced.Sincetheirintroductionin1993[23],particlefiltering techniqueshavebecomea verypopular classofnumericalmeth- odsforthesolutionofoptimalestimationproblemsinnon-linear non-Gaussian scenarios. Recently, particle based sampling filters have been proposed and used successfully to recursively update the posterior distribution using sequential importance sampling andresampling[18,19].

TheparticlefilterisanSMCmethodusedforBayesianfiltering.

Particleswithcorrespondingweightsareusedtoformanapprox- imationof a probability densityfunction (PDF). The particles are propagated over time by Monte Carlo simulation to obtain new particles and weights (usually as new information are received), henceforming a seriesof PDF approximationsover time [18,19].

SMCmethods,particlefilteringisbasedonanon-linearstatespace representationandisamethodused forstate estimation.Particle filteringisatechniqueforimplementingforrecursiveBayesianfil- teringbyusingMonte Carlobased samplingtechniques.The idea

isto approximatetheposteriordensitybya setofrandom parti- cles associated withweights. After defining theseparticles using a proposal distribution the posterior density estimates are com- putedwhichare basedontheseparticlesamplesandweights.So it is also applicable and a powerful method for filtering and/or predictionoftime seriesproblems.Particle filteringhandlesnon- linear models withnon-Gaussian noise. SMC methods provide a MonteCarlo basedsamplingmethods initself. Theyapproximate thetarget probabilitydistribution(e.g.amplitude ofspeechsignal andstockmarketreturns).Thismethodisalsoknownasbootstrap filter[17],thesurvivalofthefittest[24]orthecondensationalgo- rithm[25].

The state sequenceandthe observationsequence isa Markov randomprocess.We modelthestate (transition)equation ass_k= fs(sk−¹,uk) where sk is the state vector at time instant k, fs is the state transition function and u_k is the process noise with a known distribution. The observation sequence is modeled as yk= ^fy(sk,vk)whereykrepresentstheobservationvectorattime instantk, fy istheobservationfunctionandv_k istheobservation noisewithaknowndistribution.

4.1. SMCmethodsandBayesianfiltering

The particleapproximation oftheposterior distributionisup- dated recursively by propagating and updating the particles ac- cording tofiltering andpredictive densities. In orderto calculate theestimationweusedrecursiveBayesianfilteringtechniques.We first apply filtering to calculate the Filtering Density usingBayes Rules in (8) and calculate the Prediction Density by integrating transitionprobabilityandfilteringdensitiesin(9).

In summary Particle Filtering can be seen as a two steps process. This paper first representsthe posterior distribution ex- pressedasamixturedistributionofparticless_k⁽ⁱ⁾ withcorrespond- ingweights w⁽_kⁱ⁾.Whenanewobservationoccurs asetofnewN samplesfromthedistributionandthenewobservationisusedto evaluatethenewweights w⁽_kⁱ₊⁾₁usingp(yk+¹|^s(_kⁱ₊⁾₁).Afterthisstep aresamplingstepisincludedwheretheparticlesareresampledin ordertoduplicatetheparticleswithhighweightsandremovethe particleswithlowweights. Thisstepenablestheparticlefiltering approachtoremovedegeneracyoftheparticleapproximation.

4.2. SequentialImportanceSampling(SIS)&Resampling(SIR)

SequentialImportance Sampling isthe basis ofthe Sequential Monte Carlo(SMC) methods. The idea is to first update filtering densityusingBayesianfiltering andthen compute integralsusing importancesampling.

Hereinthispaperweusedthe priorastheproposaldistribu- tion s_k⁽ⁱ⁾∼^p(s_k|^sk−¹) and weight computation and normalization are doneasfollows:

ˆ

w⁽_kⁱ⁾

=

w⁽_kⁱ₋⁾₁p

(

yk

|

s⁽_kⁱ⁾

)

(19)

w⁽_kⁱ⁾

=

^w

ˆ

⁽_kⁱ⁾



_i^N₌₁^w

ˆ

⁽_k^{i) (20)}

But, the SIS algorithm has some limitations. The variance of theimportanceweightsincreasestochasticallyovertime[26].This variance increase posesproblemsbecausethe proposaldensityis preferred to be as close as to the posterior density and causes weight degeneracy with harmful effects on the accuracy of the simulations. Calculating the estimation on the same weight and thesample set alsowastesthecomputationresources. Toaddress thisdegeneracyproblemofSIS,aselection,resamplingprocess is usedtoeliminateparticleswithlowweightsandmultiplythepar- ticles with highimportant weights. There are different sampling

approachesincludingresidual,systematic, multinomialresampling andstratifiedsampling. InSIR,particles s_k⁽ⁱ⁾are replicatedinpro- portiontotheirweightsw⁽_kⁱ⁾andproduceN newsamplesallwith equalweights1/N,{^s(_kⁱ⁾,w⁽_kⁱ⁾}^Ni=¹.

4.3. AuxiliaryParticleFiltering(APF)

SIR approach oftenwastes a lot ofcomputation time ongen- erating particles that have small weights. In order to solve this weightdegeneracyproblem,wecanalternativelyincorporateother informationwhichwill putmoremasson relevantparticles.This should improve how our algorithm performs. This is the moti- vation for auxiliary particle filters (APF), developed by Pitt and Sheppard(1999).

Instead of applying a blind approach, we aim to update the mostpromisingparticles.Theidea hereisto usesomeformpre- dictivepowersuchasthelikelihoodinform p(y_k+1|

μ

_k+1)inpar- ticleselectionprocess toremove theweightdegeneracyproblem.

We setup a

μ

⁽_kⁱ₊⁾₁ asan estimate ofthe mode ofnewstate s_k+1, the mode of the transition density p(s_k+1|^s_k⁽ⁱ⁾) for each particle.

i=¹,2,3,. . .N. The quality of the estimate is evaluated by the auxiliaryweights:

δ

⁽_kⁱ₊⁾₁

=

^w(_k^i)p

(

y_k+1

| μ

⁽_kⁱ₊⁾₁

)

(21) Alarge auxiliary weight indicates betterrepresentation power forthe relatedparticles in accordanceto the underlyingprocess.

Forparticle selection, auxiliary weights are normalizedand they are replaced by new samples with probabilities proportional to δ_k⁽ⁱ₊⁾₁ where{^j}^{isthenewindicesand}{^N}^{representsthesizeof} thenewparticleset.Wekeptthesize ofthesample setconstant {^N}^{asthesizeoftheparticlesetsize} {^N} ^{throughoutthispaper} andwilluse N as thenewparticle setsize aswell butitisalso possibleto changethesize ofthesample set ateachtime point.

Thegenerated setofindices {^j}^{represents thesetofmostlikely} pathsofthevolatilityprocesswilltakegiventhearrivalofthenew observationy_k+1.Wewillusethenewselectedsamplesettogen- eratethenewparticless⁽_k+^j)₁∼^p(sk+¹|^s(_k^j)).New auxiliaryweights areevaluatedasfollows:

w_k+1

=

^p

(

y_k+1

|

s⁽_k+^j)₁

)

p

(

y_k+1

| μ

_k^(j₊⁾₁

)

(22)

5. MixtureKernelSmoothing(MKS)basedparameterfiltering

Statefiltering problemis turnedinto ageneraljointstate and parameter filtering problem. In this paper, we applied Liu and West’sMixtureKernelSmoothingapproachforapproximatingthe densities[Liu–West2001].We willfirstintroduce theregime pa- rameter values at time k, θ_k⁽ⁱ⁾, and provide a joint importance sample set of state and regime parameter values for the model {^s(_kⁱ⁾,θ_k⁽ⁱ⁾}i^N=¹ with associated weights {^w(_kⁱ⁾}i^N=¹. The aim is con- verted to approximate a jointposterior p(s_k+1,θ|^y1:^k) instead of justaposteriorforthestates.Wecandecomposethejointposte- riorintofollowingthreefactors:

p

(

s_k+1

, θ |

^y1:k+1

) =

^p

(

y_k+1

|θ,

^s1:k+1

)

p

(

s_k+1

|θ,

^y1:k

)

p

(θ |

^y1:k

)

(23) Thisis the product of the marginal likelihoodgiven the state andparameters,predictionoflatentvariablesgiventheparameters andpastdataandtheposteriordensityfortheregimeparameters given the past data respectively. If we know the parameter val- ues,Eq.(24)simplifiesjusttothestatefilteringproblem.Herewe

follow Liu–West (2001) and implement Mixture Kernel Smooth- ing(MKS)approachforapproximatingthedensityofp(θ|^y1:^k).We will haveMonteCarlosamples,regime parameters θ_k⁽ⁱ⁾ andasso- ciatedweights{

ω

⁽_kⁱ⁾}^N_i=1 ^{attimek fromtheposterior}distributions.

Weapproximatetheparameterdensitybyasmoothkerneldensity asfollows:

p

(θ |

y_1:k

) = 

_i^N₌₁w_k⁽ⁱ⁾N

(θ |

m⁽_kⁱ⁾

,

b²V_k

)

(24) where N(θ|^m,V)isamultivariate normaldistributionwithmean vector m andcovariancematrix V .Kernelsmoothingmethodap- proximates theposter densityfor θ by amixture of multivariate normal distributions andthe weights that comes fromthe parti- cles.b² isthekernelshrinkageparameterinLiuandWest(2001) andisusuallychosentobe (0,1)andusually decreasesslowlyas the particlesetsize increasessothat allthe parameterestimates areconcentratedclosertothemean.V_kistheparticlesetvariance andθkistheparticlesetmeanandtheyarecomputedasfollows:

V_k

= 

_i^N₌₁^w(_kⁱ⁾

(θ

_k⁽ⁱ⁾

− ¯θ

k

)(θ

_k⁽ⁱ⁾

− ¯θ

k

)

^T (25)

¯θ

_k

= 

_i^N₌₁^w_k⁽ⁱ⁾

θ

_k⁽ⁱ⁾ (26)

LiuandWest’sshrinkagemethodisusedtoremovedegeneracy problemwhichwilloccur ifm⁽_kⁱ⁾ isconsidered tobejustθ_k⁽ⁱ⁾.Liu and West took m⁽_kⁱ⁾=^aθ_k⁽ⁱ⁾+ (¹−^a) ¯θk where a=√

1−^{b2. This} mean value makes the centres of Gaussian centres to be closer andasaresultthedistributionwillhavethintails.HeandMaheu (2010)usedadiscountfactorδ∈ (⁰,1)tocontroltheshrinkageof thekernelmeansasb²=¹− [³δ−¹/2δ]^2anda=√

1−^{b2 [10].}

6. AuxiliaryParticleFilteringbasedMultipleChangepoint DetectioninGARCH/EGARCHmodels

Inthispaper,APFiscombinedwithkernelsmoothingtodesign a SMC approach for MCD problems in financial time series. The proposed methodisanAPFbasedBayesianinference approachto theMCDinfinancialtimeseriesusingvolatilitymodelsofGARCH andEGARCH.

The parameters θ are re-parameterized to take values on the real line since they will be sampled through a mixture of nor- mal kernels. The reparameterization approach is explained in the experiments sections in details. We will use θ for the re- parameterized values. The number of regimes and changepoints intheseriesisnotknownsincewe aredealingwiththeproblem of estimatingGARCH and EGARCHbased volatility models which aresubjecttoanunknownnumberofchangepointsusingSequen- tial MonteCarlo methods.A maximumnumber, an upper bound for thenumber ofchangepoints is definedwhich alsorepresents thenumberofstates,M.¯

We studied over severalalternatives includingpartial change- point models whereonly apartof theparameters are subjectto change or full changepoint models where all the parameters are subjectto change.We used both GARCH andexponentialGARCH models as the volatility processes of the time series in thispa- per.TheEGARCHmodelallowsgoodnewsandbadnewstohavea differentimpact onvolatilityonthe otherhandGARCH doesnot.

Badnews hasagreater impact onthe volatilitythan goodnews.

TheEGARCHmodelalsoallowsbignewstohaveagreaterimpact on volatility than GARCH model. We first used Normal distribu- tion inourstudiesbutempiricalstudiesoftensuggest fat tailsin the distributionof asset returns,therefore,we extended thepro- posedsolutiontoincorporatestudent-tdistributionwith

ν

degrees offreedomaswell.Wedefineallthesealternativemodelsandthe algorithmproposedbelow.

Fig. 7. States & filtered changepoints, log returns and volatility respectively.

6.1.Auxiliaryparticlefilteringbasedmultiplechangepointdetection algorithm

The proposed approach is used to approximate p(s_k,θ|^y1:^k) givenasetofparticles{^s_k⁽ⁱ⁾,θ_k⁽ⁱ⁾}^N_i=1^andweights{^w_k⁽ⁱ⁾}_i^N₌₁ ^foreach regimewhereN isthesizeoftheparticleset.Latentvariables sk, andparameterestimationsaredoneusingtheSMCapproachinthe algorithmbelow.Latentvolatilities

σ

_k² areinferredthroughGARCH andEGARCHrecursionsusingEqs.(12)and(14).

Thealgorithmoftheproposedapproachisgivenbelow:

For∀^i,i=¹,2,. . .N wedefine

μ

⁽_kⁱ₊⁾₁ asthemodeofthepredic- tionofthelatentvariable p(sk+¹|^s(_kⁱ⁾,θ_k⁽ⁱ⁾).

1. Produce stratified uniform random variables {^ui}ⁿ_i=1 ^within uniformintervalsu_i∼

(ⁱ⁻_n¹,_nⁱ)

2. Compute the auxiliary weights δ_k⁽₊ⁱ⁾₁ ∼

ω

⁽_kⁱ⁾p(y_k+1|^y₁⁽ⁱ_:⁾_k,m⁽_kⁱ⁾) where m⁽_kⁱ⁾=^aθ_k⁽ⁱ⁾+ (¹−^a) ¯θ_k and ¯θ_k⁼_i^N₌₁

ω

_k⁽ⁱ⁾θ_k⁽ⁱ⁾ and draw a sample fromi=¹,2,. . .N with thecorresponding auxiliary weights.

3. Sample the newregime parameter set from the kth normal componentofthekernelfunctiongiven:

θ

_k⁽ⁱ₊⁾₁

∼

^N

(

m_k⁽ⁱ⁾

,

b²V_k

)

Vk

= 

_i^N₌₁w_k⁽ⁱ⁾

(θ

_k⁽ⁱ⁾

− ¯θ

k

)(θ

_k⁽ⁱ⁾

− ¯θ

k

)

^T b²

=

1

− [

3

δ −

1

/

2

δ]

²

4. Sampleanewstatevalues_k⁽ⁱ₊⁾₁∼^p(s⁽_k+ⁱ⁾₁|^s(_kⁱ⁾,θ_k⁽₊ⁱ⁾₁) 5. Correcttheweights

ω

_k⁽₊ⁱ⁾_1∼ ^p(yk+1|s

(i) k+1,θ_k⁽ⁱ⁾₊₁) p(y_k+1|μ⁽ⁱ⁾_k+1,m⁽ⁱ⁾_k)

6. Repeat from 2nd to 5th steps using stratified sampling ap- proachdefinedbelow.Stratified samplinghereisusedtosta- bilize estimates over multiple runs [10]. Stratified Sampling is based on ideas used in survey sampling consists in pre- partitioning the (0,1] ^{interval into n disjoint sets. To pro-} duce a new sample size of {ⁿ} ^{from a population} {^sk}_k^K₌₁ withweights{^wk}_k^K₌₁^{,stratifiedsamplingfirstproducesstrati-} fieduniformrandomvariables{^ui}ⁿ_i=1^{withinuniformintervals} u_i∼

(ⁱ⁻_n¹,_nⁱ).u_i’saredrawnindependentlyineachofthese sub-intervals.

7. Experiments&results

Toanalyzetheperformanceoftheproposedapproach,wesetup twosimulationexperiments,onewithsyntheticdataandtheother one withan emerging markettime series,BIST100which isde- scribed inSection 2, usingGARCH andEGARCH volatilitymodels withboth Gaussian andStudent-t distributions and usingpartial andfullparameterchangepointmodelswhereselectedparameters oralltheparameters ofthevolatilitymodels arechangingindif- ferentregimes.WealsocomparedtheresultsoftheproposedAPF basedmultiplechangepointdetectionalgorithmwithHe andMa- heu’s GARCH based model and reported results below [10]. The

Fig. 8. Top left: log returns, top right: changepoint estimates, center left: true volatility, center right: volatility estimates, bottom left: effective sample size (ESS).

resultsforthesyntheticdataandrealtime seriesdataarerepre- sentedrespectively.

7.1. Experimentsonsyntheticdataseries

Wegeneratedthesynthetictime seriesdatausinganEGARCH modelwithGaussian distributionandfull parameterchangepoint modeling.This time seriesincludesswitching regimesby provid- ing different

ω

k,

α

k, βk and

γ

k parameters at changepoint loca- tions. Synthetictime seriesis generated usinga predefined left–

rightHMM model.The parameter configuration ofthe algorithm includesthemaximumnumberofthechangepoints,initialparam- eters for each regime and volatility model parameter sets. First the maximum limit for the number of changepoints (the num- ber ofregimes)is defined. Thishasnot to bethe numberofthe exact number of changepoints but will be defined asthe maxi- mum numberof changepoint limit that the algorithmcan reach.

Fourdifferentregimeparametersaredefinedforeachregime,θm=

[

ω

k,

α

k,βk,

γ

k]^{. The upperbound forthe numberof statesis de-} finedasM=^{5.Statetransition}probabilitiesaredefinedasP .The priorparametersforeachregimearerandomlygeneratedas

ω

k∼ Gamma(1,0.5),

α

kwithBeta(4,1),β_kwithBeta(1,8)∗ (¹−

α

k)and

γ

k withBeta(1,8)foreachregime.GammaandBetadistributions areusedtohaveparametersbetween(0,1).θm isreparameterized aslog(θm).

Thisresultswithprobabilitiescloseto1whichmeansahigher probability for not changing the states and not producing too manyregimes.Atthetimewherethetransitionprobability islow enough to move to next state, we create a changepoint location andcreatethechangepointtime serieswitha500time step.We define regimes as the series between two changepoint locations andfor eachregime thevolatility modelparameters are selected to be low if it isodd, high ifthe regime numberis even in or- der to see thefluctuations inthe synthetictime series.We then calculatethevolatility usingEGARCHmodelwithGaussiandistri- bution. Changepointlocations, volatilityandthe generatedreturn

Fig. 9. Filtered parameters in each regime (ωk,αk, βk,γk) and their %2,5 and %97,5 quantiles.

seriesare plotted inFig. 7 below. The blue dashed vertical lines state the exact changepoint locations which allow us to easily identifythestateswitchingregimes.Wetried severalalternatives regardingtheparameterstoseetheperformanceofthealgorithm.

In thesimulated data set several deltaparameters are evaluated such as 0.25, 0.5, 0.75, 0.9 and 0.99 and selected as the best to be used in this paper. We had experiments with a 500time stepsyntheticdataserieswithdifferentparticlesizesof500,1000 and 5000.

Fig. 7presentstherealchangepointlocationsforfivedifferent statesandregimesinblueandcorrespondingfilteredchangepoints and estimated regimes in red. The second and third subfigures showsthe generated log returns of thetime series andEGARCH basedvolatilityrespectivelywithbluedashedverticallinesstating thetruechangepointlocationsinthecorrespondingseries.

We reported the results and estimates of the proposed APF basedMCDalgorithminFig. 8.Thetopleftandrightfiguresgive the real logreturn seriesand corresponding filtered changepoint locations using the volatility model.Here in the top rightfigure we will easily identifythat the algorithmis able toestimate the changepointlocationsverycloselytorealregimeswitchingpoints.

The middle part presents the filtered volatility estimates along withthetruevolatilities. We alsoplotthesequence oftheeffec- tivesamplesize(ESS)ofparticlestocheckweightdegeneracy.ESS ofparticlesisequaltoNwhenallweightsareequalto1/N.Inthe plotthereexistsalargefluctuationaroundregimeswitchingpoints andoutliersandtheyreturntonormalquicklyaftertheoutlieror theregimeswitchingpoint.

Fig. 9presentsthefilteredparametersandtheir%2.5 and%97.5 quantiles inred andgreen and their mean values as blacklines