[R] Estimation of AR(1) Model with Markov Switching
napps22
n.j.apps22 at gmail.com
Thu Dec 1 17:58:00 CET 2011
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.
More information about the R-help
mailing list