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. ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.