In mhunter1/dynr: Dynamic Models with Regime-Switching

Introduction

knitr::opts_chunk$set(echo = TRUE,message=FALSE,warning=FALSE)

This example illustrates how to fit a discrete-time regime-switching linear model in the dynr package.

Data

First, create a dynr data object using dynr.data. Here we have 6 observed variables.

require(dynr)

data(NonlinearDFAsim)
data <- dynr.data(NonlinearDFAsim, id="id", time="time",observed=colnames(NonlinearDFAsim)[c(3:8)])

Measurement Model

Next, we specify the measurement model using prep.measurement. The first three of the six observed variables load on the Positive Emotion (?) latent variable, and the last three load on the Negative Emotion variable. Parameters are indicated by parameter names (e.g.lambda_* in the following code), and fixed values are indicated by "fixed". The values.* arguments specify the starting values and fixed values. The params.* arguments specify the free parameter names or the reserved word "fixed" for fixed parameters. Parameters with the same name are constrained to be equal. Since we do not have any covariate in this model. No *.exo arguments are supplied. Dispite this being a regime-switching model, we assume the same measurement model in the two regimes. Hence only one measurement model needs to be specified.

$\begin{bmatrix}{} y1(t) \ y2(t) \ y3(t) \ y4(t) \ y5(t) \ y6(t) \ \end{bmatrix} = \begin{bmatrix}{} 1 & 0 \ \lambda_{21} & 0 \ \lambda_{31} & 0 \ 0 & 1 \ 0 & \lambda_{52} \ 0 & \lambda_{62} \ \end{bmatrix} \begin{bmatrix}{} PE(t) \ NE(t) \ \end{bmatrix} + \epsilon$,$\epsilon\sim N\Big( \begin{bmatrix}{} 0 \ 0 \ 0 \ 0 \ 0 \ 0 \ \end{bmatrix} , \begin{bmatrix}{} \epsilon_{1} & 0 & 0 & 0 & 0 & 0 \ 0 & \epsilon_{2} & 0 & 0 & 0 & 0 \ 0 & 0 & \epsilon_{3} & 0 & 0 & 0 \ 0 & 0 & 0 & \epsilon_{4} & 0 & 0 \ 0 & 0 & 0 & 0 & \epsilon_{5} & 0 \ 0 & 0 & 0 & 0 & 0 & \epsilon_{6} \ \end{bmatrix}$

meas <- prep.measurement(
  values.load=matrix(c(1, .8, .8, rep(0, 3),
                       rep(0, 3), 1, .8, .8), ncol=2),
  params.load=matrix(c("fixed", "lambda_21", "lambda_31", rep("fixed",3),
                       rep("fixed",3), "fixed", "lambda_52","lambda_62"), ncol=2),
  state.names=c('PE', 'NE'))

Regime Probabilities

In the next step, we specify transition probability matrix between the regimes using prep.regimes. This transition probabilities from time t to time t+1 are:

| | Regime1(t+1) | Regime2(t+1) | | ------------ |:-------------:|:------------:| | Regime1(t) | $$ \frac{exp(p11)}{exp(p11)+exp(0)} $$ | $$ \frac{exp(0)}{exp(p11)+exp(0)} $$ | | Regime2(t) | $$ \frac{exp(p21)}{exp(p21)+exp(0)} $$ | $$ \frac{exp(0)}{exp(p21)+exp(0)} $$ |

(Here, $p11$ and $p21$ are model parameters. They can be perceived as odd ratios. )

regimes <- prep.regimes(
    values=matrix(c(.9, 0, 0, .9), 2, 2), 
    params=matrix(c("p11", "fixed", "fixed", "p22"), 2, 2))

Dynamic Model

In the next chuck, we specify our dynamic models by first specifying the covariance matrices of measurement errors and dynamic noises using prep.noise, and then specifying the dynamic functions using prep.formulaDynamics. The dynamic models are:

\begin{align} \text{Regime 1:}&\ &PE(t+1) = a1 \times PE(t) + w1(t),\ &NE(t+1) = a2 \times NE(t) + w2(t),\ &w(t) \sim N\Big( \begin{bmatrix}{} 0.00 \ 0.00 \ \end{bmatrix} , \begin{bmatrix}{} \zeta_{1} & 0 \ 0 & \zeta_{2} \ \end{bmatrix} \Big)\ \text{Regime 2:}&\ &PE(t+1) = a1 \times PE(t) + c12 \times \frac{\exp(abs(NE(t)))}{1 + \exp(abs(NE(t)))} \times NE(t) + w1(t),\ &NE(t+1) = a2 \times NE(t) + c21 \times \frac{\exp(abs(PE(t)))}{1 + \exp(abs(PE(t)))} \times PE(t) + w2(t),\ &w(t) \sim N\Big( % Fri Sep 23 13:24:01 2016 \begin{bmatrix}{} 0.00 \ 0.00 \ \end{bmatrix} , \begin{bmatrix}{} \zeta_{1} & 0 \ 0 & \zeta_{2} \ \end{bmatrix} \Big) \end{align}

We assume the same dynamic noise process applies to both regimes, hence we only have one matrix for dynamic noise specification (*.latent) in prep.noise. The *.observed arguments in prep.noise specify the measurement error ($\epsilon$ in the measurement model above).

mdcov <- prep.noise(
    values.latent=diag(0.3, 2),
    params.latent=diag(paste0("zeta_",1:2), 2),
    values.observed=diag(0.1, 6),
    params.observed=diag(paste0("epsilon_",1:6), 6))

We write dynamic functions in a list that contains two lists of functions, one for each regime.

formula=list(
  list(PE~a1*PE,
       NE~a2*NE),
  list(PE~a1*PE+c12*(exp(abs(NE)))/(1+exp(abs(NE)))*NE,
       NE~a2*NE+c21*(exp(abs(PE)))/(1+exp(abs(PE)))*PE) 
)

Since this is a nonlinear model we need to specify the jacobian matrix containing the differentiation of each dynamic function above with respect to each latent variable (PE \& NE).

jacob=list(
  list(PE~PE~a1,
       NE~NE~a2),
  list(PE~PE~a1,
       PE~NE~c12*(exp(abs(NE))/(exp(abs(NE))+1)+NE*sign(NE)*exp(abs(NE))/(1+exp(abs(NE))^2)),
       NE~NE~a2,
       NE~PE~c21*(exp(abs(PE))/(exp(abs(PE))+1)+PE*sign(PE)*exp(abs(PE))/(1+exp(abs(PE))^2))))

Then we combine dynamic functions and the jocobian matrix together with parameter specifications in prep.formulaDynamics.

dynm<-prep.formulaDynamics(formula=formula,startval=c(a1=.3,a2=.4,c12=-.5,c21=-.5),isContinuousTime=FALSE,jacobian=jacob)

The prep.tfun function is used here to transform regime switching probability parameters (odd ratios) to transition probabilities. This function can also be used to tranform parameters on a constrained scale to an unconstrained scale (e.g. exponential transformation to ensure parameters take positive values).

trans<-prep.tfun(formula.trans=list(p11~exp(p11)/(1+exp(p11)), p22~exp(p22)/(1+exp(p22))), formula.inv=list(p11~log(p11/(1-p11)),p22~log(p22/(1-p22))), transCcode=FALSE)

Initial Values

After that, we specify values at time t=0 using prep.initial. These values are used to initialize the recursive algorithm (extended Kalman filter) that dynr uses. The *.inistate arguments specify the initial (starting) states of the latent state variables. The *.inicov arguments specify the starting covariance matrix of the latent state variables. The *.regimep specifies the initial probabilities of the two regimes.

initial <- prep.initial(
    values.inistate=c(0, 0),
    params.inistate=c("fixed", "fixed"),
    values.inicov=diag(1, 2), 
    params.inicov=diag("fixed", 2),
    values.regimep=c(.8, .2),
    params.regimep=c("fixed", "fixed")
)

Dynr Model

Now we put together everything we've previously specified in dynr.model. This code connects the recipes we've written up with our data and writes a c file in our working directory. We can inspect c functions that go with each recipe in the c file.

model <- dynr.model(dynamics=dynm, measurement=meas, noise=mdcov, 
                    initial=initial, regimes=regimes, transform=trans, 
                    data=data, 
                    outfile="RSNonlinearDiscrete")

Tex Options

We can check our model specifications in a neatly printed pdf file using the following code.

The printex command is used to write the model into a Latex file, with a name given by the outFile argument. Then, the tools::texi2pdf command generates a pdf file from the latex file we just created. The system command prints out the pdf file:

We can also print out the model in R, instead of generating a Latex file, using the command plotFormula.

printex(model,ParameterAs=model$param.names,printInit=TRUE, printRS=TRUE,
        outFile="RSNonlinearDiscrete.tex")
tools::texi2pdf("RSNonlinearDiscrete.tex")
system(paste(getOption("pdfviewer"), "RSNonlinearDiscrete.pdf"))