[R] Estimation of AR(1) Model with Markov Switching

[napps22]

Dear R users,

I have been trying to obtain the MLE of the following model

state 0: y_t = 2 + 0.5 * y_{t-1} + e_t
state 1: y_t = 0.5 + 0.9 * y_{t-1} + e_t

where e_t ~ iidN(0,1)

transition probability between states is 0.2

I've generated some fake data and tried to estimate the parameters using the
constrOptim() function but I can't get sensible answers using it. I've tried
using nlminb and maxLik but they haven't helped. Any tips on how I could
possibly rewrite my likelihood in a better way to improve my results would
be welcome.

Given below is my R code

# markov switching regime model
# generate data for a AR(1) markov switching model with the following pars
# state 0: y_t = 2 + 0.5 * y_{t-1} + e_t
# state 1: y_t = 0.5 + 0.9 * y_{t-1} + e_t
# where e_t ~ N(0,1) 
# transition probabilities p_s0_s1 = p_s1_s0 = 0.20

# generate realisations of the state

gamma_s0 <- qnorm(0.8)
gamma_s1 <- qnorm(0.2)
gamma <- rep(0,100)
state <- rep(0,100)

# choose initial state at t=0 to be state 0 

gamma[1] <- gamma_s0
state[1] <- 0

for(i in 2:100) {
	
	if(rnorm(1) < gamma[i-1]) {
	gamma[i] <- gamma_s0
	state[i] <- 0
	}
	else {
	gamma[i] <- gamma_s1
	state[i] <- 1
	}
	}
	
# generate observations
# choose y_0 = 0
# recall state at t=1 was set to 0 

y1 <- 2 + 0.5 * 0 + rnorm(1)
y <- rep(0,100)
y[1] <- y1

for(i in 2:100) {
	
	if(state[i]==0) {
		y[i] <- 2 + 0.5 * y[i-1] + rnorm(1)
		}
	else {
		y[i] <- 0.5 + 0.9 * y[i-1] + rnorm(1)
		}
		}

# convert into time series object
 
y <- ts(y, start = 1, freq = 1)

# construct negative conditional likelihood function 

neg.logl <- function(theta, data) {

# construct parameters	
    beta_s0 <- theta[1:2]
	beta_s1 <- theta[3:4]
	sigma2 <- exp(theta[5])
	gamma0 <- theta[6]
	gamma1 <- theta[7]
	
# construct probabilities

#probit specification
  
p_s0_s0 <- pnorm(gamma_s0)
p_s0_s1 <- pnorm(gamma_s1)
p_s1_s0 <- 1-pnorm(gamma_s0)
p_s1_s1 <- 1-pnorm(gamma_s1)

# create data matrix

X <- cbind(1,y)

# assume erogodicity of the markov chain
# use unconditional probabilities
 
p0_s0 <- (1 - p_s1_s1) / (2 -p_s0_s0 -p_s1_s1)
p0_s1 <- 1-p0_s0

# create variables
 
p_s0_t_1 <- rep(0, nrow(X))
p_s1_t_1 <- rep(0, nrow(X))
p_s0_t <- rep(0, nrow(X))
p_s1_t <- rep(0, nrow(X))
f_s0 <- rep(0,nrow(X)-1)
f_s1 <- rep(0,nrow(X)-1)
f <- rep(0,nrow(X)-1)
logf <- rep(0, nrow(X)-1)

p_s0_t[1] <- p0_s0
p_s1_t[1] <- p0_s1 

# initiate hamilton filter

for(i in 2:nrow(X)) {
	
# calculate prior probabilities using the TPT
# TPT for this example gives us 
# p_si_t_1 = p_si_t_1_si * p_si_t + p_si_t_1_sj * p_si_t
# where p_si_t_1 is the prob state_t = i given information @ time t-1
# p_si_t_1_sj is the prob state_t = i given state_t_1 = j, and all info @
time t-1
# p_si_t is the prob state_t = i given information @ time t

# in this simple example p_si_t_1_sj = p_si_sj

p_s0_t_1[i] <- (p_s0_s0 * p_s0_t[i-1]) + (p_s0_s1 * p_s1_t[i-1])
p_s1_t_1[i] <- (p_s1_s0 * p_s0_t[i-1]) + (p_s1_s1 * p_s1_t[i-1])

# calculate density function for observation i
# f_si is density conditional on state = i
# f is the density 
  
f_s0[i] <- dnorm(y[i]-X[i-1,]%*%beta_s0, sd = sqrt(sigma2)) 
f_s1[i] <- dnorm(y[i]-X[i-1,]%*%beta_s1, sd = sqrt(sigma2))
f[i] <- (f_s0[i] * p_s0_t_1[i]) + (f_s1[i] * p_s1_t_1[i])

# calculate filtered/posterior probabilities using bayes rule
# p_si_t is the prob that state = i given information @ time t

p_s0_t[i] <- (f_s0[i] * p_s0_t_1[i]) / f[i]
p_s1_t[i] <- (f_s1[i] * p_s1_t_1[i]) / f[i]

logf[i] <- log(f[i])

}

logl <-sum(logf)
return(-logl)

}


# restrict intercept in state model 0 to be greater than intercept in state
model 1
# thus matrix of restrictions R is [1 0 -1 0 0 0 0]

R <- matrix(c(1,0,-1,0,0,0,0), nrow = 1)

# pick start values for the 7 unknown parameters

start_val <- matrix(runif(7), nrow = 7)

# ensures starting values are in the feasible set 

start_val[1,] <- start_val[3,] + 0.1

# estimate pars 

results <-constrOptim(start_val,neg.logl,grad = NULL, ui = R, ci = 0)

Regards,

N

--
View this message in context: http://r.789695.n4.nabble.com/Estimation-of-AR-1-Model-with-Markov-Switching-tp4129417p4129417.html
Sent from the R help mailing list archive at Nabble.com.