Markov-Switching Regressions Part 1 – Manual R Code for the basic quasi-maximum likelihood (univariate case)

I am sorry that I haven’t posted any new blog entries for the last couple months. I have been busy with my work and my side project which I do not want to discuss openly. Anyway,  I would like to talk about Markov-switching regression today.

When we run time series regressions,  it is usually the case where our estimated coefficients  may not be stable and be subject to regime shifts. (This can be checked by CUSUM squared tests and other break point tests).  To model this kind of regime shift , several types of econometric models were developed by Econometricians. One of them is the Markov-switching type regressions, which has gained recognition ever since the influential work by Hamilton (1989) . Nowadays, we have many statistical packages to run markov switching univariate or multivariate regressions. However,  I never found anyone trying to post the manual code for the markov-switching regressions. So, let us write the code together! 🙂

The code I am going to present is the Quasi-likelihood Estimation method for a dynamic 2-state markov-switchhing AR(1) model:

Let Y_t be the observed variables. Assume our model has the following form:

 Y_t = \alpha_s+ \beta_s Y_{t-1}+\epsilon_t

where s= 0,1 signals 2 different states. \alpha_s and \beta_s are subject to regime shift depending on the state changes

First, Assuming asymptotic normality, and let \Theta = (\alpha_s, \beta_s), we can write the density of Y_t conditional on Y_{t-1} and s_t=i  as:

f(Y_t|s_t, Y_{t-1}; \Theta)= \frac{1}{\sqrt{2\pi\sigma_\epsilon^2}}exp{\frac{-\epsilon_{st}^2}{2\sigma_\epsilon}}

Second, we have: f(Y_t|Y_{t-1}; \Theta)=P(s_t = 0|Y_{t-1};\Theta) f(Y_t|s_t=0,Y_{t-1};\Theta)+P(s_t=1|Y_{t-1};\Theta)f(Y_t|s_t=1,Y_{t-1};\Theta)

Third, we have: P(s_t=i|Y_t; \Theta)=\frac{P(s_t=i|Y_{t-1};\Theta)f(Y_t|s_t=i; \Theta)}{f(Y_t|Y_{t-1}; \Theta)}

The fourth equation we have is: P(s_{t+1}=i| Y_t; \Theta) = p_{0i}P(s_t=0|Y_t; \Theta)+p_{1i}P(s_t=i|Y_t;\Theta)

where p_{0i} and p_{1i} are the transition probabilities from previous state to the next state. For further details of the Markov-switching model, please refer to most papers by Hamilton (1989,1994, 2005) etc.. Here, I only briefly outlined the key equations to set up the likelihood equation

Now if we are given initial values of P(s_t = 0|Y_{t-1};\Theta), we iterate the 4 equations to obtain the sum of the likelihood. Then the quasi-likelihood equation is:

L(\Theta)= \frac{1}{T}\sum_{t=1}^{T}{lnf(y_t|Y_{t-1};\Theta)}

To maximize the equation above, we can use the optimx package in R or the nlm function in R

Here I post my R code for setting up the quasi-likelihood function for the example above:

likeli_fun <-function(param){

y <- data[,1] #Here, use your own dataset, dependent variable, put it in column 1
alpha<-param[1] #Intercept in regime 0
alpha1 <- param[2] #intercept in regime 1
var <-param[3]^2 # variance of the error
p00 <- param[4] #The transition probability
p11 <-param[5]
beta <-param[6] #AR(1) slope in regime 0
beta1 <-param[7] #AR(1) slope in regime 1

#Set initial value for prob1 and prob0:
A <- matrix(c(1-p00,-1+p11,1,-1+p00,1-p11,1),nrow=3,ncol=2) #This is the A matrix mentioned in Hamilton 1994

inv_A <-solve(t(A)%*%A)
sol_A <- inv_A%*%t(A)  #Why we do this? Please check out Hamilton (1994)

prob0<-sol_A[1,3] #Recommended initial values by Hamilton (1994)

#Start iterating to get the likelihood
T <- length(y)
likv <-0
j_iter <-2
while (j_iter <= N) {

eps0 <- y[j_iter]-alpha-beta*y[j_iter-1]
eps1 <- y[j_iter]-alpha1-beta1*y[j_iter-1]
f0 <- (1/sqrt(2*pi*var))*exp(-0.5*eps0*eps0/var) #equation 1 for regime0
f1 <- (1/sqrt(2*pi*var))*exp(-0.5*eps1*eps1/var) #equation 1 for regime1

con_fz <- prob0*f0+prob1*f1  #Equation 2

fil_prob0 <- (prob0*f0)/con_fz   #Equation 3
fil_prob1 <- (prob1*f1)/con_fz

prob0 <- p00*fil_prob0+(1-p00)*fil_prob1 #Equation 4
prob1 <- (1-p11)*fil_prob0+p11*fil_prob1
lik <- log(con_fz)   # Take the log
likv <- likv-lik     # add to the last likelihood
j_iter <- j_iter+1   #Start another round of iteration

Next step, all you need to do in R is to fire up nlm or optimx in R to maximize the likelihood to compute the estimates of the parameters.

In my next post, let us start using the packages in R and command in Eviews to do some interesting analysis using Markov-switching regressions.


Cointegration Tests that allow for structural breaks?

Structural breaks can sometimes be a huge hassle to deal with when you are trying to investigate the long run relationship between the variables by running standard cointegration tests such as Engle-Granger or the standard Johansen test.

Now, when the structural breaks are present, one solution is the residual-based cointegration test proposed by Gregory and Hansen (1996).   The relevant R programs and example can be found on website of Bruce E. Hansen by clicking here. Note that p is number of variables and r is number of cointegrating rank being tested

The test proposed by Johansen et al. (2000) appears to be another solution. The relevant R program for computing the critical values can be found at Dave Giles’ website here. (In the program, remember to set the correct breakpoint proportion and the value of q!!)

To do the Johansen et al. (2000) test, it can be decomposed into the following steps:

Step1: Identify the structural breaks.

Step2:  incorporate the date dummy (D2), trend*dummy interaction term(D2*@trend), as well as the shift indicator dummy(I2)  as exogenous variables into the original VAR. In Dave Giles example,  the variables he adds are  D2(-2), I2(-1), trend*D2(-2), and I2.

step 3: construct the usual Johansen trace statistics

(How to calculate the trace statistics? See this paper here )

The asymptotic critical values depend on the proportion of the way through the sample that the break occurs (λ = 0.44 in our case); and on (pr), where p is the number of variables under test p = 2, here), and r is the cointegrating rank being tested. So, for us,  r  = 0, 1.  Unlike Gregory-Hansen (1996) test, the Johansen et al. (2000) test can be modified to allow for two structural breaks.


I downloaded the EViews workfile,  which is available here


Now Time for some sleep before tomorrow’s work….Zzzzzzz