# DLM package in R to estimate a state space model with drifts

I am trying to use the DLM package in R to estimate a state space model where the measurement and transition equations are as follows.

The measurement equations are: \begin{align} \left( \begin{array}{c} p_t \\ n_t \end{array} \right) = \left( \begin{array}{c} \psi p_{t-1} \\ \phi n_{t-1} \end{array} \right) + \left( \begin{array}{cc} 1 & -\psi \\ 1-\phi & 0 \end{array} \right) \left( \begin{array}{c} v_t \\ v_{t-1} \end{array} \right) + \left( \begin{array}{c} \epsilon_t \\ \omega_t \end{array} \right) \end{align}

and the transition equations are: \begin{align} \left( \begin{array}{c} v_t \\ v_{t-1} \end{array} \right) = \left( \begin{array}{cc} 1 & 0 \\ 1 & 0 \end{array} \right) \left( \begin{array}{c} v_{t-1} \\ v_{t-2} \end{array} \right) + \left( \begin{array}{c} r_t \\ 0 \end{array} \right) \end{align}

with $$\epsilon_t, \omega_t, r_t$$ being independent random errors with non-zero means.

I have two questions. One is fundamental to state space models and the other is about the coding in R.

1. How to model the non-zero means? Are these included in the parameter matrices or in the state vector? Given that they are constant it seems to be that the former is preferred.

2. Why does this code not work? Following the above choice about modelling the means, my likelihood function (using dlm) looks like this.

###########################################################################
# DLM PACKAGE

my_dlmfc <- function(par=c(psi,phi,mu_p,mu_n,mu_v,sig_p,sig_n,sig_v)){

psi = par[1]
phi = par[2]

mu_p = par[3]
mu_n = par[4]
mu_v = par[5]

sig_p = par[6]
sig_n = par[7]
sig_v = par[8]

# Z is the state vector and Y is the vector of the observations.
# Note, I have expanded the state vector to include $$p_{t-1}$$ and $$n_{t-1}$$, and the scalar 1:

#(#1) Z(t) = (p(t-1) n(t-1) v(t) v(t-1) 1)'
#(#2) Y(t) = (p(t) n(t))'

# The measurement equations are
#(#3) Z(t) = GG Z(t-1) + w(t)

# The tramsition equations are
#(#4) Y(t) = FF Z(t) + v(t)

GG <- matrix(c(0, 0, 0, 0,    0,
0, 0, 0, 0,    0,
0, 0, 1, 0, mu_v,
0, 0, 1, 0,    0,
0, 0, 0, 0,    1),
nrow=5,byrow=TRUE)

W  <- diag(c(0, 0, sig_v^2, 0, 0))

FF <- matrix(c(psi,   0,     1, -psi, mu_p,
0, phi, 1-phi,    0, mu_n),
nrow=2,byrow=TRUE)

V  <- diag(c(sig_p^2,sig_n^2))

## starting values
m0 <- matrix(c(4, 4, 4, 4,1))
C0 <- diag(c(10,10,10,10,0))

my_dlm <- dlm(FF=FF, V=V, GG=GG, W=W, m0=m0, C0=C0)

return(my_dlm)
}


My MLE estimation using dlmMLE then is:


#par=c(psi,phi,mu_p,mu_n,mu_v,sig_p,sig_n,sig_v)){
mleold <- dlmMLE(obs,
c(0.8, 0.8, 5, 5, 5, 1, 1, 1),  #set initial values for parameters
my_dlmfc, method = "BFGS")

comp <- data.frame(row.names=c("psi","phi","mu_p","mu_n","mu_v","sig_p","sig_n","sig_v","convergence"), MLE = c(mleold$$par,mleold$$convergence))
knitr::kable(comp,digits=4, caption = "MLE estimates")


The results look "reasonable" to me:

Table: MLE estimates

|            |    MLE|
|:-----------|------:|
|psi         | 0.7742|
|phi         | 0.7643|
|mu_p        | 4.8983|
|mu_n        | 4.9967|
|mu_v        | 4.9626|
|sig_p       | 1.0040|
|sig_n       | 0.8518|
|sig_v       | 1.0527|
|convergence | 0.0000|


I then go on an estmiated the filtered values of the state vector $$v_t$$ using: