I did not quite get what the problem was from your description ...
However, I did search on CRAN and found at least two packages that can fit
markov switching ar models so that may be an easier way to go.
Hth, Ingmar
On Thu, Dec 1, 2011 at 5:58 PM, napps22 <n.j.apps22@gmail.com> wrote:
> 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.
>
[[alternative HTML version deleted]]