This example illustrates how to fit a discrete-time regime-switching linear model in the dynr package.
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)])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'))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", 0, 0, "p22"), 2, 2))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)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")
)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")We can check our model specifications in a neatly printed pdf file using the following code. The printex() function produces a latex file with model specifications. The ParameterAs argument
printex(model,ParameterAs=model$param.names,printInit=TRUE, printRS=TRUE,
outFile="RSNonlinearDiscrete.tex")
tools::texi2pdf("RSNonlinearDiscrete.tex")
system(paste(getOption("pdfviewer"), "RSNonlinearDiscrete.pdf"))Finally, it is time to cook dynr (i.e. fit our model through parameter optimization)!
res <- dynr.cook(model)And serve!
summary(res)## names parameters s.e. t-value ci.lower
## a1 a1 0.21984699 0.019003498 11.568764 0.18260082
## a2 a2 0.25822127 0.019660186 13.134223 0.21968801
## c12 c12 -0.68483308 0.082124861 -8.338925 -0.84579485
## c21 c21 -0.64524920 0.086870577 -7.427707 -0.81551240
## lambda_21 lambda_21 1.19031882 0.023219942 51.262782 1.14480857
## lambda_31 lambda_31 1.19883399 0.023239888 51.585189 1.15328464
## lambda_52 lambda_52 1.09744529 0.019643179 55.869028 1.05894537
## lambda_62 lambda_62 0.93392064 0.017156953 54.433945 0.90029363
## zeta_1 zeta_1 0.33646588 0.014703861 22.882825 0.30764684
## zeta_2 zeta_2 0.27983488 0.010729134 26.081776 0.25880616
## epsilon_1 epsilon_1 0.27943138 0.008362887 33.413267 0.26304042
## epsilon_2 epsilon_2 0.10685304 0.006009070 17.781957 0.09507547
## epsilon_3 epsilon_3 0.09513918 0.005940007 16.016679 0.08349698
## epsilon_4 epsilon_4 0.12040683 0.004794840 25.111751 0.11100911
## epsilon_5 epsilon_5 0.12134955 0.005345098 22.702960 0.11087335
## epsilon_6 epsilon_6 0.10544462 0.004191294 25.158008 0.09722983
## p11 p11 0.96327332 0.013024302 73.959687 0.93774616
## p22 p22 0.83713991 0.046808018 17.884541 0.74539788
## ci.upper
## a1 0.2570932
## a2 0.2967545
## c12 -0.5238713
## c21 -0.4749860
## lambda_21 1.2358291
## lambda_31 1.2443833
## lambda_52 1.1359452
## lambda_62 0.9675477
## zeta_1 0.3652849
## zeta_2 0.3008636
## epsilon_1 0.2958223
## epsilon_2 0.1186306
## epsilon_3 0.1067814
## epsilon_4 0.1298045
## epsilon_5 0.1318258
## epsilon_6 0.1136594
## p11 0.9888005
## p22 0.9288819
##
## -2 log-likelihood value at convergence = 28362.81
## AIC = 28398.81
## BIC = 28506.92