# Specifying a State Space Model with the MARSS package

I would like to use the state space technique outlined in the paper by Mokinskia et al. (https://www.american.edu/cas/economics/research/upload/2013-4.pdf) to estimate time varying and asymmetric threshold values in the context of the Carlson-Parkin approach to quantify ordered variables. The description of the concept can be found on pages 4 and 5, the specification of the model to be estimated can be found on pages 6 to 8. I tried various packages like dlm, KFAS but eventually went with the "MARSS" package as it seems to be the only one which allows me to specify the model exactly as it is outlined in the paper (the upper part at page 8) and uses the same Kalman Filter implementation. The results are however strange as you will see later. The following setup was used:

The data set required to run the code below (in data.table ):

test = structure(list(ue = c(7.8, 7.7, 7.5, 7.5, 7.5, 7.2, 7.5, 7.4,
7.4, 7.2, 7.5, 7.5, 7.2, 7.4, 7.6, 7.9, 8.3, 8.5, 8.6, 8.9, 9,
9.3, 9.4, 9.6, 9.8, 9.8, 10.1, 10.4, 10.8, 10.8, 10.4, 10.4,
10.3, 10.2, 10.1, 10.1),
vdo = c(0.816, 1.148, 1.581, 3.54, -3.319,
1.608, 1.247, 1.344, 1.201, 1.192, 1.434, 1.377, 1.489, 1.41,
1.213, 0.893, 0.952, 0.768, 0.863, 0.996, 0.926, 1.106, 1.071,
1.136, 1.025, 1.018, 1.346, 1.537, 2.07, 1.726, 3.026, -2.489,
-0.381, -0.451, -0.009, -1.8),
vup = c(0.184, -0.148, -0.581,-2.54, 4.319, -0.608, -0.247, -0.344, -0.201, -0.192, -0.434,
-0.377, -0.489, -0.41, -0.213, 0.107, 0.048, 0.232, 0.137, 0.004,
0.074, -0.106, -0.071, -0.136, -0.025, -0.018, -0.346, -0.537,
-1.07, -0.726, -2.026, 3.489, 1.381, 1.451, 1.009, 2.8)),
.Names = c("ue","vdo", "vup"), row.names = c(NA, -36L), class = c("data.table","data.frame"))


The variables ue, vdo and vup are $\mu$, as well as the upper and the lower row of the LHS matrix in equation 5 in the paper.

require(MARSS)
require(data.table)
p          = 1
m          = 2
TT         = nrow(test)
y          = matrix(as.matrix(test[,list(ue)]),nrow=p)  #y
Z          = array(t(test[,list(vdo,vup)]),c(p,m,TT))   #data matrix - x in the paper
U          = "zero"                                     #intercepts
A          = "zero"                                     #intercepts
B          =  array(list(1,0,0,1),c(m,m,1))             #T in paper
Q          =  array(list("q",0,0,"q"),c(m,m,1))         #Q in paper ("diagonal and equal")
R          =  array(list("r"),c(p,p,1))                 #H in paper
A          =  "zero"
x0         =  array(c(1,2),c(m,1,1))                    #initial values
V0         =  "diagonal and unequal"                    #initial values
model.list = list(Z=Z,U=U,A=A,B=B,Q=Q,R=R,x0=x0,V0=V0,tinitx=0)
fit$par = fit$start
res        = cbind(matrix(kfl$xtT,ncol=2),t(y))  The model runs but the following observations make me think if I did something wrong: The estimated state vectors are along my understanding supposed to be the upper and the lower thresholds around$\mu$. However in some cases (as you can see when you run the code and look at the result matrix res)$\mu\$ is outside of the estimated interval. Moreover when you use the whole data set (I just supplied three years of data here for the purpose of the convenience of readers) there are also periods in which the lower threshold is above the upper threshold which is obviously wrong. What am I am doing wrong? I am glad for all hints! Thank you very much!