A regime switching time series model of daily river flows Krisztina - - PowerPoint PPT Presentation

A regime switching time series model of daily river flows

Krisztina Vasas 26.08.2005. 14th European Young Statisticians Meeting, Debrecen

Characteristics of river flow series

– Hydrologists do not like simulated river flows, because they often do not feature

• the highly skewed marginal distribution of the empirical series
• easily separable ascending and descending periods (regimes)
• the asymmetric shape of the hydrograph (i.e. that short

ascending periods are followed by long descending ones)

– Moreover, in empirical river flow series:

• The ascending and descending regimes have random durations
• The discharges have random increments
• The dynamics are different in the two distinct regimes

Why do we like regime switching models?

• A linear model – even if generated with skewed and

seasonal driving force – is not appropriate for approximating – the high quantiles – the probability density – the assymmetric shape of the empirical river flow series

• A possible alternative is to use a light-tailed conditionally

heteroscedastic model – While it solves the first two problems, the assymetry of regime lengths remains unresolved

• A natural approach: regime switching models
• Obvious regime specification: periods with increasing and

decreasing discharges

The model of Lu and Berliner (1998)

• Three states (they call them: normal, rising and falling)
• Every state has an autoregressive structure with normally

distributed noise

• The lagged precipitation is included as an additional

regressor in the rising regime

• We do not have precipitation data, and if precipitation is

not included, the distribution of this model is symmetric

• A possible way out: non-symmetric generating noise in

the autoregressive regime switching model

A regime switching model

• It indicates the actual regime-type:

– 0 if t belongs to the ascending regime – 1 if t belongs to the descending one.

• It is a Markov-chain of hidden states with the transition

matrix

• The generating noises are conditionally independent valued:

– ε1,t ~ Γ(α,λ) in the ascending regime (shocks to the system) – ε2,t ~ N(0, σ2) in the descending regime

( )

   = + + − = + =

− −

1 , ,

, 2 1 , 1 1 t t t t t t t

I if c c Y a I if Y Y ε ε

        =

11 10 01 00

p p p p P

Stationarity of the model

• The model can be written as a stochastic difference

equation:

• where
• It follows from Brandt (1986) that a unique stationary

solution is :

• as

t t t t

b Y a Y + =

−1

{ } { }

1

1

= = +

=

t t

I I t

a a χ χ

{ } { }

1 , 2 , 1 = = +

=

t t

I t I t t

b χ ε χ ε

∑

∞ = − − − −

+ =

1 1 1 i i t i t t t t t

b a a a b Y L

( ) ( )

( ) ( )

∞ < <

+

log and , log b E a E

Model estimation

• Parameters of the model: θ = (p00, p11,α, λ, η, a, c),

where η=1/σ2 .

• The latent variables (It) make model estimation

complicated.

• Obvious choices:

– Maximum likelihood (The likelihood can be written down, but only recursively because of the latent variables) – Markov Chain Monte Carlo (MCMC) – Efficient Method of Moments (EMM).

Posterior distribution in the Bayesian framework

{ } { }

( )

{ } { }

( ) { } ( )

= ∝

= = =

) ( | , | | ,

1 1 1

θ θ θ θ f I f I Y f Y I f

T t t T t t t t T t t

( )

[ ] {

}[

] {

}

) ( )) ( (

11 10 01 00 2

11 10 01 00 1 1 1 ) , ( 1 ) , (

θ

χ σ χ λ α

f p p p p c Y a c Y f Y Y f

n n n n T t I t t N I t t

t t

⋅ − − − − =∏

= = − = − Γ

{ } ( ) { } { }

( ) { }

T t t t T t t t

I d Y I f Y f

1 1

| , |

= =

∫

= θ θ

where for example n00=#(It-1=0 and It=0).

Markov Chain Monte Carlo method

• We need to sample from the joint posterior density to get

the estimation of the parameters

• The method for this is to create a Markov-chain with

stationary distribution .

• The hidden states (It) is updated as well as the structural

parameters

• Gibbs-sampling is used as much as possible
• By the help of conjugate priors the full conditionals can be

computed for six out of the seven parameters and for the hidden states

• A Metropolis-Hastings step is necessary for the shape

parameter (α) of the increments in the ascending regime { } { }

( )

t T t t

Y I f | ,

1 θ =

Metropolis-Hastings

• Sample a candidate point Y from a proposal distribution

q(y|θ(t)) at each iteration t

• Accept the candidate point with probability

then θ(t+1)=y and any other case θ(t+1)= θ(t).

( )

( )

( )

( )

( )

( )

( )

( )

( ) 

       = 1 , | | min ,

t t t t

y q y q y y θ θ π θ π θ α

Gibbs sampling

• If you can write down the full conditionals:
• the sampling can be realized in that way:

( ) ( ) ( )

∫

+ − + − + −

=

j n j j j n j j j n j j j

dθ θ θ θ θ θ π θ θ θ θ θ π θ θ θ θ θ π , , , , , , , , , , , , , , , , , |

1 1 1 1 1 1 1 1 1

K K K K K K

( )

( )

( )

( )

( )

( )

) 1 ( 1 ) 1 ( 2 1) (t 1 1 ) ( ) ( 3 1) (t 1 2 1 2 ) ( ) ( 3 ) ( 2 1 1 1

, , , θ | ~ , , , θ | ~ , , , | ~

+ − + + + + + + t n t n t n t n t t t n t t t

θ θ θ π θ θ θ θ π θ θ θ θ θ π θ K M K K

Choice of priors

• α ~ Γ(αu,λu)
• λ~ Γ(r,β)
• η ~ Γ(q,ρ)
• a ~ N(µ,τ)
• c~ N(ν,κ)
• p00 ~ β(u1,v1)
• p11 ~ β(u2,v2)

Full conditionals for λ and η

• distribution of λ:
• distribution of η:

{ }{ } ( )

) ) ( , ( ~ , , |

1

β α α λ + − + Γ

∑

= −

t

I t t t t

Y Y r n I Y

{ } { } ( ) ( ) ( )

        + − − − + Γ

∑

= − 1 2 1 1

, 2 ~ , , , |

t

I t t t t

c Y a c Y q n c a I Y ρ η

Full conditionals for a and c

• distribution of c:
• distribution of a:

{ } { } ( ) ( ) ( ) ( ) ( )

            + − + − + − − ∑

= −

κ κ υκ

2 1 2 1 1 2 1

1 1 , 1 1 ~ , , | a n a n aY Y a N a I Y c

t

I t t t t

{ } { } ( ) ( )( ) ( ) ( ) 

           − + − + − − +

∑ ∑ ∑

= − = − = − 1 2 1 1 2 1 1 1

1 , ~ , , |

t t t

I t I t I t t t t

c Y c Y c Y c Y N c I Y a η τ η τ η µτ

Full conditionals for the transition probabilities

• distribution of p00:
• distribution of p11:

{ } ( ) ( )

1 01 1 00 00

, ~ | v n u n I p

t

+ + β

{ } ( ) ( )

1 10 1 11 11

, ~ | v n u n I p

t

+ + β

Full conditionals for the hidden states

( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( )( ) ( ) ( )

) ( 1 1 ) ( , , 1 , 1 | 1 ) ( 1 1 , , , |

1 2 11 1 11 00 1 2 11 1 1 1 1 11 00 1 2 00 1 2 00 1 1 1

c Y a c Y N p Y Y p p c Y a c Y N p Y Y I I I P c Y a c Y N p p Y Y p Y Y p Y Y I I I P

t t t t t t t t t t t t t t t t t t t t t t

− − − + − Γ − − − − − = = = = − − − − − + − Γ − Γ = = = =

− − − − + − − − − − + −

Updating the parameter α

• In the Metropolis-Hastings step for α we use a normal

distribution for updating: α*~ N(α,δ2)

• Symmetry implies a simple acceptance ratio:
• where π denotes the posterior distribution

( )

( ) 

       α π α π

*

, 1 max

Results for River Tisza

• Observation period: 10 years (3650 days)
• 30000 updates for all parameters
• Only the last 15000 members of the chain will

participate in the analysis

• Various starting values and hyperparameters

were tried to get information about the stability

• f parameters

Results for α and λ

• Posterior mean of α is 1.003, indicating that the increments of

the rising regimes are close to exponential

• Average height of the increment in the rising regime:

α/λ=106.68 m3/s

Results for η

• Posterior mean of the standard deviation of the noise in the

descending regime is 25.83 m3/s

Results for a and c

• Posterior mean of a is 0.815, indicating high level of persistence

even in the descending regime

• Posterior mean of c is 104.9

Results for the transition probabilities

• Average duration (calculated from the posterior means
• f the transition probabilities)

– of the ascending period is 2.38 days – of the descending period is 14.12 days

• indicating much longer descending than ascending

regimes

Convergence properties

• All parameters apart from α stationarize after at

most 1000 iterations

• Parameter α reaches equilibrium after a longer

period (to determine its stationary distribution more accurately, more iterations might be needed)

• Acceptance ratio for α in the Metropolis-Hastings

algorithm is 0.51

A trajectory of the original and of the simulated process

Density of the original (red) and simulated (blue) series

What if the state process is not Markovian?

• As a consequence of Markovity the time spent in one state

is geometrically distributed.

• A possible generalization:

– The length of ascending periods is distributed as NegBin(b,p1) – The length of descending periods is distributed as Geom(p0) – Regime durations are independent

• Special case: b=1 (we obtain the previous model)

Estimation of the generalized model

• The structural parameters are: θ = (α, λ, η, a, c, p0, b, p1)
• Now we cannot determine conditional probability for the

hidden states

• Instead, we introduce the change points of the regimes as

latent variables:

– end of descending periods: s=(s1,...,sm) – end of ascending periods: t=(t1,...,tm)

• Starting values and number of the change points are created

from the results of the previous model

• If we know the change points, {It} can be determined, so the

algorithm remains almost the same as previously

• The posterior distribution:
• With an appropriate choice of priors six out of eight

structural parameters and the change points can be updated with Gibbs-sampling

• Two parameters need Metropolis-Hastings-steps:

– parameter α of the increment in the ascending regime – parameter b of the length of the ascending regime

• Priors are the same as in the previous model, the b

parameter has a Gamma-distributed prior

) ( ) ( ) ( ) , , | ( ) ( ) | , ( ) , , | ( ) | , , (

1 ) ( Geom 1 ) , ( NegBin

1

θ θ θ θ θ θ f t s f s t f s t Y f f s t f s t Y f Y s t f

i i p m i i i p b

        − − = = ∝

− =

∏

Updating the change points

• Full conditional for the s :
• where the second terms are : (N1 is negative binomial, N2 is

geometrically distributed)

{ } ( ) =

+ =

− − t i i i i

Y t t x t s P , , |

1 1

{ }

( ) (

) { }

( ) (

)

∑

− − − − ≤ ≤ − − − − ≤ ≤

+ = + = + = + =

− −

y i i i i i i i i t t t t i i i i i i i i t t t t

t t y t s P y t s t t Y P t t x t s P x t s t t Y P

i i i i

, | , , | , | , , |

1 1 1 1 1 1 1 1

1 1

( ) ( ) ( ) =

− = + = = + =

− − − 1 2 1 1 1 1

, |

i i i i i i

t t N N P x N P t t x t s P

( ) ( ) ( ) ( )

∑

− − − − − − − −

− − y i i i i

p y t t Geom p b y NegBin p x t t Geom p b x NegBin

11 1 00 11 1 00

1 , 1 , , 1 , 1 , ,

Updating the change points

• Full conditional for the t :
• where the second terms are : (N1 is negative binomial, N2 is

geometrically distributed)

{ } ( ) =

+ =

+ t i i i i

Y s s x s t P , , |

1

{ }

( ) (

) { }

( ) (

)

∑

+ + ≤ ≤ + + ≤ ≤

+ = + = + = + =

+ +

y i i i i i i i i s t s t i i i i i i i i s t s t

s s y s t P y s t s s Y P s s x s t P x s t s s Y P

i i i i

, | , , | , | , , |

1 1 1 1

1 1

( ) ( ) ( ) =

− = + = = + =

+ + i i i i i i

s s N N P x N P s s x s t P

1 2 1 2 1,

|

( ) ( ) ( ) ( )

∑

− − − − − − − −

+ + y i i i i

p b y s s NegBin p y Geom p b x s s NegBin p x Geom

00 1 11 00 1 11

1 , , 1 , 1 , , 1 ,

Estimation results

• Posterior mean of b is 1.1244
• b=1 cannot be rejected on a 5

% level

• Other parameters remain

largely unaffected

• To alter the number of

change points reversible jump MCMC is needed

Posterior density of b

References

• Brandt,A. (1986) : The stochastic equation Yn+1=AnYn+Bn with

stationary coefficients, Advances in Applied Probability, 18 211-220

• Elek, Péter - Márkus, László (2004): A long range dependent model

with nonlinear innovations for simulating daily river flows, Natural Hazards and Earth Systems Sciences, 4, 277-283.

• Lu, Zhan-Qian és Berliner, L.Mark (1999): Markov Switching time

series models with application to a daily runoff series, Water Resour. Res., 35(2), 523-534, 1999.

Thank you for your attention!