Regime Switching State Space Model with R code

[This article was first published on K & L Fintech Modeling, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

This post explains a Markov regime switching state space model. The bottom line is two-fold: 1) expanding states by each regime transitions and 2) collapsing each updated estimates for the next state prediction. The step 2) is necessary to fix the dimension of previous states so that Kalman recursion does not produce exponentially increasing states. This combined filter is called Kim filter (= Kalman filter + Hamilton filter + Kim collapsing procedure).


Regime Switching State Space Model



In the previous posts below, we have explained and implemented regime switching models.


In this post, to complete a series of posts regarding the regime switching model, a regime switching state space model is investigated. The content of this post is heavily borrowed from Kim (1994) and Kim and Nelson (1999).



Regime Switching State Space Model


As an example of a regime switching state space model, Prof. Kim used the following generalized Hamilton model for the log of real GNP (Lam; 1990) in his paper and book.

ln(GNPt)=nt+xtnt=nt1+μ0+μ1stxt=ϕ1xt1+ϕ2xt2+ututN(0,σ2)st=0,1Pij=[p00p01p10p11]

where ln(GNPt) is a real GNP level. nt is a deterministic series with a regime-switching growth rate and xt is stationary AR(2) cycle process.

Since ln(GNPt) is the log level variable, the difference of it, yt=ln(GNPt)ln(GNPt1), can be represented as a state space model in the following way.

yt=μ0+μ1st+xtxt1xt=ϕ1xt1+ϕ2xt2+ut

It is a typical approach that a state space model is represented as a vector-matrix form for using Kalman filter as follows.

yt=μst+[11][xtxt1][xtxt1]=[ϕ1ϕ210][xt1xt2]+[ut0]utN(0,σ2)μst=μ0+μ1st,μ1>0st=0,1Pij=[p00p01p10p11]yt=μst+Fxtxt=Axt1+vt

For the sake of notational simplicity, we use xt instead of xt.



Kalman Filtering


Kalman filter with regime switching is used to get state estimates from a state space model taking regime transition into account and has the following recursion.

xijt|t1=Axit1|t1Pijt|t1=APit1|t1AT+Qηijt|t1=ytμjFxijt|t1Hijt|t1=FPijt|t1FT+RKij=Pijt|t1FT[Hijt|t1]1xijt|t=xijt|t1+Kijηijt|t1Pijt|t=(IKijF)Pijt|t1

In regime-dependent Kalman filter, all the notations are augmented with superscript {ij} except xit1|t1 and Pit1|t1 since these two estimates are in i-state (two-state) but other estimates must reflect state transitions from i to j (four-state). For example, xt|t1 and xijt|t1 are different in terms of conditioning information.

xt|t1=E[Xt|ψt1]xijt|t1=E[Xt|ψt1,St=j,St1=i]

In contrast to the single regime, however, in the multiple regimes, xijt|t and Pijt|t cannot be used the next state prediction due simply to the mismatch both 1) between xijt|t and xit1|t1 and 2) between Pijt|t and Pit1|t1. To resolve this mismatch problem, Kim (1994) developed a dimension collapsing algorithm.



Kim(1994)’s Collapsing procedure


Kim (1994) introduces a collapsing procedure (approximation) to reduce the (M×M) posteriors (xijt|t and Pijt|t) into M to complete the above Kalman filter recursion.

xjt|t=Mi=1P[St1i,Stj|ψt]xijt|tP[St=j|ψt]Pjt|t=Mi=1P[St1i,Stj|ψt][Pijt|t+(xjt|t – xijt|t)(xjt|t – xijt|t)T]P[St=j|ψt]

To calculate the above approximation, when we calculate P[St1i,Stj|ψt], P[St=j|ψt] is easily obtained by summing its M branches from each i states.

P[St=j|ψt]=Mi=1P[St1i,Stj|ψt]

Knowing P[St1i,Stj|ψt] means that we observe time t data since the last time of information set is t and the state is migrated from i to j. For this data information into account, we can think of the marginal probability of state transition by integrating out data.

P[St1i,Stj|ψt]=f(yt,St1i,Stj|ψt1)f(yt|ψt1)f(yt|ψt1)=Mj=1Mi=1f(yt,St1i,Stj|ψt1)

As can be seen from the above equations, when we know f(yt,St1i,Stj|ψt1) which is the joint density of data and two states, P[St1i,Stj|ψt] and f(yt|ψt1) are easily obtained. We get the following the joint density.

f(yt,St1i,Stj|ψt1)=f(yt|St1i,Stj,ψt1)×P(St1i,Stj|ψt1)

Now we need to know two parts. The first part is the forecast error given data. (MVN : probability density function of multivariate normal distribution with zero mean and forecast error variance)

f(yt|St1i,Stj,ψt1)=MVN(forecast error,its variance)

The second part is calculated by the multiplication of the transition probability and the summation of its branches.

P(St1i,Stj|ψt1)=P[Stj|St1i]×Mk=1f(St2k,St1i|ψt1)=Pij×f(St1=i|ψt1)

In the above equation, as we know, Pij is already known as transition probability matrix and f(St2k,St1i|ψt1) is exactly what we want to find but evaluated at the previous t-1 time. Therefore for the iteration, f(S1k,S0i|ψ0) calls for initialization with the steady state probabilities.

Therefore, we can calculate P[St1i,Stj|ψt] and P[Stj|ψt] through the above equations.



Kim (1994) Filter for Regime Switching State Space model


Kim filtering procedure is summarized in a sequence of equations.

Kalman Filtering
xijt|t1=Axit1|t1Pijt|t1=APit1|t1AT+Qηijt|t1=ytμjFxijt|t1Hijt|t1=FPijt|t1FT+RKij=Pijt|t1FT[Hijt|t1]1xijt|t=xijt|t1+Kijηijt|t1Pijt|t=(IKijF)Pijt|t1
Hamilton Filtering
f(yt,St1i,Stj|ψt1)=N(ηijt|t1,Hijt|t1)×Pij×P(St1i|ψt1)f(yt|ψt1)=Mj=1Mi=1f(yt,St1i,Stj|ψt1)P[St1i,Stj|ψt]=f(yt,St1i,Stj|ψt1)f(yt|ψt1)P[St=j|ψt]=Mi=1P[St1i,Stj|ψt]
Kim’s Collapsing
xjt|t=Mi=1P[St1i,Stj|ψt]xijt|tP[St=j|ψt]Pjt|t=Mi=1P[St1i,Stj|ψt][Pijt|t+(xjt|t – xijt|t)(xjt|t – xijt|t)T]P[St=j|ψt]


In particular, three red colored terms (xit1|t1, Pit1|t1, and P(St1i|ψt1)) are initialized for iteration to be started. Once the iterations get started, these red colored terms are replaced with blue colored terms (xjt|t, Pjt|t, and P(Stj|ψt)) for each iterations.



R code


The following R code estimate parameters of example model in Kim (1994) with U.S. real GNP data ranging from 1952.4 to 1984.4.

Regime Switching State Space Model Kim filter R code

The original gauss code can be found at Prof. Kim’s web site (http://econ.korea.ac.kr/~cjkim/SSMARKOV.htm) where you can find bunch of gauss programs for regime switching state space model and Gibbs sampling. In the following R code, I use the same variable name in most case for a comparison purpose.



#========================================================#
# Quantitative ALM, Financial Econometrics & Derivatives 
# ML/DL using R, Python, Tensorflow by Sang-Heon Lee 
#
# https://kiandlee.blogspot.com
#——————————————————–#
# Lam’s generalized Hamilton model in Kim (1994)
#========================================================#
 
graphics.off()  # clear all graphs
rm(list = ls()) # remove all files from your workspace
 
#=======================================
#             Functions
#=======================================
 
# parameter constraints
trans < function(b0){
    b1 < b0
    
    # probability
    b1[1:2< exp(b0[1:2])/(1+exp(b0[1:2]))
    
    # variance
    b1[5< b0[5]^2
    
    # AR(1), AR(2) coefficients
    XX6 < b0[6]/(1+abs(b0[6]))
    XX7 < b0[7]/(1+abs(b0[7]))
    b1[6< XX6 + XX7
    b1[7< 1*XX6*XX7
    
    return(b1)
}
 
# Kim filter = Hamilton + Kalman
loglike < function(param_unc, y) {
  
    # param_unc = init_param_unc; y = y
    nobs < length(y)
    
    param < trans(param_unc)
    
    p11 < param[1# Pr[St=1/St-1=1]
    p00 < param[2# Pr[St=0/St-1=0]
    MU0 < param[3]       # delta0
    MU1 < MU0 + param[4# delta0 + delta1
    VAR_W < param[5]     # sigma^2
    phi1 < param[6]      # AR(1)
    phi2 < param[7]      # AR(2)
 
    R < matrix(0,1,1)
    Q < rbind(c(VAR_W, 0), c(00))
    H < t(c(11))
    F < rbind(c(phi1, phi2), c(10))
    
    # initialization
    B_LL0 < B_LL1 < param[8:9#c(0,0)
    P_LL1 < P_LL0 < 
      matrix(solve(diag(4 kronecker(F,F))%*%
               as.vector(Q), nrow = 2)
    
    # initialization with steady state probabilities
    PROB_1 < (1p00)/(2p11p00) #Pr[S_0=1/Y_0]
    PROB_0 < 1PROB_1            #Pr[S_0=0/Y_0]
      
    #——————————————
    #  START ITERATION
    #——————————————
    likv < 0.0
    for(t in 1:nobs) {
      
        #———————————
        # KALMAN FILTER
        #——————————— 
            
        # state prediction
        B_TL00 < F %*% B_LL0  # S=0, S’=0
        B_TL01 < F %*% B_LL0  # S=0, S’=1
        B_TL10 < F %*% B_LL1  # S=1, S’=0
        B_TL11 < F %*% B_LL1  # S=1, S’=1
                            
        # state prediction uncertainty
        P_TL00 < F %*% P_LL0 %*% t(F) + Q
        P_TL01 < F %*% P_LL0 %*% t(F) + Q
        P_TL10 < F %*% P_LL1 %*% t(F) + Q
        P_TL11 < F %*% P_LL1 %*% t(F) + Q
        
        # forecast   
        fit00 < H %*% B_TL00 + MU0
        fit01 < H %*% B_TL01 + MU1
        fit10 < H %*% B_TL10 + MU0
        fit11 < H %*% B_TL11 + MU1
        
        # forecast error
        F_CAST00 < y[t]  fit00
        F_CAST01 < y[t]  fit01
        F_CAST10 < y[t]  fit10
        F_CAST11 < y[t]  fit11
                    
        # forecast error variance
        SS00 <  H %*% P_TL00 %*% t(H) + R
        SS01 <  H %*% P_TL01 %*% t(H) + R
        SS10 <  H %*% P_TL10 %*% t(H) + R
        SS11 <  H %*% P_TL11 %*% t(H) + R
        
        # Kalman gain
        kg00 < P_TL00 %*% t(H) %*% solve(SS00)
        kg01 < P_TL01 %*% t(H) %*% solve(SS01)
        kg10 < P_TL10 %*% t(H) %*% solve(SS10)
        kg11 < P_TL11 %*% t(H) %*% solve(SS11)
        
        # updating
        B_TT00 < B_TL00 + kg00 %*% F_CAST00
        B_TT01 < B_TL01 + kg01 %*% F_CAST01
        B_TT10 < B_TL10 + kg10 %*% F_CAST10
        B_TT11 < B_TL11 + kg11 %*% F_CAST11
                    
        P_TT00 < (diag(2 kg00 %*% H ) %*% P_TL00
        P_TT01 < (diag(2 kg01 %*% H ) %*% P_TL01
        P_TT10 < (diag(2 kg10 %*% H ) %*% P_TL10
        P_TT11 < (diag(2 kg11 %*% H ) %*% P_TL11
                    
        #——————————— 
        # HAMILTON FILTER
        #——————————— 
                        
        # Pr[St,Yt/Yt-1]
        NORM_PDF00 < dnorm(F_CAST00,0,sqrt(SS00))
        NORM_PDF01 < dnorm(F_CAST01,0,sqrt(SS01))
        NORM_PDF10 < dnorm(F_CAST10,0,sqrt(SS10))
        NORM_PDF11 < dnorm(F_CAST11,0,sqrt(SS11))
        
        PR_VL00 < NORM_PDF00*(  p00)*PROB_0
        PR_VL01 < NORM_PDF01*(1p00)*PROB_0
        PR_VL10 < NORM_PDF10*(1p11)*PROB_1
        PR_VL11 < NORM_PDF11*(  p11)*PROB_1
    
        # f(y_t|Y_{t-1})
        PR_VAL < PR_VL00 + PR_VL01 + PR_VL10 + PR_VL11
        
        # Pr[St,St-1/Yt]
        PRO_00 < as.numeric(PR_VL00/PR_VAL)
        PRO_01 < as.numeric(PR_VL01/PR_VAL)
        PRO_10 < as.numeric(PR_VL10/PR_VAL)
        PRO_11 < as.numeric(PR_VL11/PR_VAL)
                                  
        # Pr[St/Yt]
        PROB_0 < PRO_00 + PRO_10
        PROB_1 < PRO_01 + PRO_11
                                            
        #——————————— 
        # Kim’s Collapsing 
        #——————————— 
        B_LL0 < (PRO_00*B_TT00 + PRO_10*B_TT10)/PROB_0
        B_LL1 < (PRO_01*B_TT01 + PRO_11*B_TT11)/PROB_1
        
        temp00 < (B_LL0B_TT00) %*% t(B_LL0B_TT00)
        temp10 < (B_LL0B_TT10) %*% t(B_LL0B_TT10)
        temp01 < (B_LL1B_TT01) %*% t(B_LL1B_TT01)
        temp11 < (B_LL1B_TT11) %*% t(B_LL1B_TT11)
                   
        P_LL0 < (PRO_00*(P_TT00 + temp00) +
                  PRO_10*(P_TT10 + temp10))/PROB_0
        P_LL1 < (PRO_01*(P_TT01 + temp01) +
                  PRO_11*(P_TT11 + temp11))/PROB_1
                                            
        likv = likv log(PR_VAL)
    } 
    return(likv)
}
 
#=======================================
#                  Run
#=======================================
 
# Lam’s Real GNP Data Set : #1952.3 ~ 1984.4
rgnp < c(  1378.21406.81431.41444.91438.2
    1426.61406.81401.21418  , 1438.81469.6
    1485.71505.51518.71515.71522.61523.7
    1540.61553.31552.41561.51537.31506.1
    1514.21550  , 1586.71606.41637  , 1629.5
    1643.41671.61666.81668.41654.11671.3
    1692.11716.31754.91777.91796.41813.1
    1810.11834.61860  , 1892.51906.11948.7
    1965.41985.21993.72036.92066.42099.3
    2147.62190.12195.82218.32229.22241.8
    2255.22287.72300.62327.32366.92385.3
    2383  , 2416.52419.82433.22423.52408.6
    2406.52435.82413.82478.62478.42491.1
    2491  , 2545.62595.12622.12671.32734  , 
    2741  , 2738.32762.82747.42755.22719.3
    2695.42642.72669.62714.92752.72804.4
    2816.92828.62856.82896  , 2942.73001.8
    2994.13020.53115.93142.63181.63181.7
    3178.73207.43201.33233.43157  , 3159.1
    3199.23261.13250.23264.63219  , 3170.4
    3179.93154.53159.33186.63258.33306.4
    3365.13444.73487.13507.43520.4
 
# # GNP growth rate : 1952.4 ~ 1984.4
< diff(log(rgnp),1)*100
 
# initial guess
init_param_unc < c(3.5,0.0,0.4,1.7,0.5,1,0.7,1,1)
 
# optimization
mle < optim(init_param_unc, loglike, method=c(“BFGS”),
             control=list(maxit=50000,trace=2), y=y)
mle < optim(mle$par, loglike, method=c(“Nelder-Mead”),
             control=list(maxit=50000,trace=2), y=y)
 
# compare estimation results with Kim (1994) paper
est_param < cbind(c(trans(mle$par), mle$value),
    c(0.9540.4561.4572.4210.773
      1.2460.3675.2240.535176.33))
 
est_param[5,1< sqrt(est_param[5,1]) # var -> std
 
colnames(est_param) < c(“our est”“Kim est”)
rownames(est_param) < 
  c(“p00”,“p11”,“delta0”“delta1”“sigma”
    “phi1”“phi2”,“x0”,“x-1”,“likv”)
 
est_param
 
cs


The following figure shows that our estimated parameters and evaluated log likelihood are very similar to results of Kim (1994).
Regime Switching State Space Model Kim filter R code, Kim (1994) Dimension Collapsing


Concluding Remarks


This post dealt with the regime switching state space model. It seems that this regime switching modeling approach is widely and actively used in trading practice. I tried to give intuitive and sequential explanations and to implement R code for the regime switching state space model by using Kim (1994) algorithm. For more information about this modeling techniques, I refer the reader to Kim and Nelson (1999).

Kim, Chang-Jin (1994) Dynamic Linear Models with Markov-switching, Journal of Econometrics 60-1, pp. 1-22.
Kim, C.-J. and C. R. Nelson (1999) State-Space Models with Regime Switching: Classical and Gibbs-Sampling Approaches with Applications. MIT Press. (http://econ.korea.ac.kr/~cjkim/SSMARKOV.htm)


To leave a comment for the author, please follow the link and comment on their blog: K & L Fintech Modeling.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)