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.
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.
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:
finalmodel <- my_dlmfc(mleold$par)
filtered <- dlmFilter(obs, finalmodel)
and here is where I am getting confused. Here are the first 10 rows of the state vector
X1 X2 X3 X4 X5
1 4.0000000 4.000000e+00 4.000000 4.000000 1.0000000
2 -2.5274746 2.321655e+00 8.962633 4.000000 0.1864689
3 -0.9907256 1.309535e-16 10.058909 8.042542 0.3608859
4 -1.5903605 -4.828191e-17 11.935773 10.002706 0.3810106
5 -1.7481068 -3.217985e-17 13.651207 11.936853 0.3606711
6 -1.8276048 -2.635410e-17 15.081485 13.555340 0.3322520
Why does the state variable X5 change (wrongly I think) from 1.00 given the transition equation matrix being
> GG(finalmodel)
[,1] [,2] [,3] [,4] [,5]
[1,] 0 0 0 0 0.000000
[2,] 0 0 0 0 0.000000
[3,] 0 0 1 0 4.962633
[4,] 0 0 1 0 0.000000
[5,] 0 0 0 0 1.000000
with a 1 in [5,5] and the covariance matrix having correctly a zero variance for the entry 1 in the state vector
> W(finalmodel)
[,1] [,2] [,3] [,4] [,5]
[1,] 0 0 0.000000 0 0
[2,] 0 0 0.000000 0 0
[3,] 0 0 1.108264 0 0
[4,] 0 0 0.000000 0 0
[5,] 0 0 0.000000 0 0
Any insight on question #1 or #2 is highly appreciated! Thanks!!
I am adding the details of the underlying model and the citation to it. The model is taken from Madhavan and Sobczyk "Price Dynamics and Liquidity of Exchange-Traded Funds" in the Journal Of Investment Management, Vol. 14, No. 2, (2016), pp. 1–17.
It relates an Exchange Traded Fund's observed price ($p$) and observed NAV ($n$) to the unobserved, true value($v$). All values are in natural logs.
The model assumes that the unobserved, true value $v$ evolves as a randomg walk with non-zero drift as in $$ v_t = v_{t-1} + r_t $$ where $r_t \sim (\mu_r, \sigma_r^2)$.
The observed values price ($p$) and observed NAV ($n$) are related to the unobserved, true value as follows: $$ p_t = \psi p_{t-1} + v_t - \psi v_{t-1} + \epsilon_t $$ where $\epsilon_t \sim (\mu_\epsilon, \sigma_\epsilon^2)$, and $$ n_t = \phi n_{t-1} + (1 - \phi) v_t + \omega_t $$ where $\omega_t \sim (\mu_\omega, \sigma_\omega^2)$. As previously mentioned, $\epsilon, \omega, r$ are independent random errors.
After the various questions I received below, let me explain further. What follows will show that the filtered value of the state vector ought to contain a ``1" in the last position for all $t$, and hence why I am confused with the output from dlmFiltered().
To simplify matters, let me write wlg a simpler model using the notation from dlm().
Let the observed variable be $y_t = x_t$ and the unobserved variable a random walk with a non-zero constant $x_t = x_{t-1} + \mu_x + w_t$, where $w_t$ is a random error term with mean zero and variance $\sigma_w^2$.
The model can then be written in state space notation as follows. The measurement equation is
$$
y_t = F \theta_t =
\left(\begin{array}{ccc} 1 & 0 & 0 \end{array}\right)
\left(\begin{array}{c} x_t \\ x_{t-1} \\ 1 \end{array}\right)
$$
and the transition equations reflecting the random walk with non-zero constant can be written as
$$ \theta_t = \left(\begin{array}{c} x_t \\ x_{t-1} \\ 1 \end{array}\right) = G \theta_{t-1} + \left(\begin{array}{c} w_t \\ 0 \\ 0 \end{array}\right) = \left(\begin{array}{ccc} 1 & 0 & \mu_x \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{array}\right) \left(\begin{array}{c} x_{t-1} \\ x_{t-2} \\ 1 \end{array}\right) + \left(\begin{array}{c} w_t \\ 0 \\ 0 \end{array}\right) $$
Following Hamilton "Time Series Analysis" 1994, we can obtain the filtered values of $\theta_{t|t}$ using the procedure described on the pages 377-381. Note, to remove clutter, I neglect the notation $\hat{}$ to identify an estimate. We start the recursion with
$$ \theta_{1|0} = \left(\begin{array}{c} 0 \\ 0 \\ 1 \end{array}\right) $$
and (using Hamilton's notation)
$$ P_{1|0} = \left(\begin{array}{ccc} v & 0 & 0 \\ 0 & v & 0 \\ 0 & 0 & 0 \end{array}\right) $$
then following Hamilton and noting that the measurement equation has no error term, we have for the filtered value of $\theta$
$$ \theta_{1|1} = \theta_{1|0} + P_{1|0} F' \left( F P_{1|0} F'\right)^{-1}\left(y_1 - F \theta_{1|0}\right) $$
and plugging in we have
$$ \theta_{1|1} = \left(\begin{array}{c} 0 \\ 0 \\ 1 \end{array}\right) + \left(\begin{array}{c} v \\ 0 \\ 0 \end{array}\right) v^{-1} y_1 = \left(\begin{array}{c} y_1 \\ 0 \\ 1 \end{array}\right) $$
Then the updated inference for $\theta_{2|1}$ is
$$ \theta_{2|1} = G \theta_{1|1} = \left(\begin{array}{ccc} 1 & 0 & \mu_x \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{array}\right) \left(\begin{array}{c} y_1 \\ 0 \\ 1 \end{array}\right) = \left(\begin{array}{c} y_1 + \mu_x \\ y_1 \\ 1 \end{array}\right) $$
and the udpated MSE $P_{2|1}$ is given by
$$ P_{2|1} = G [P_{1|0} - P_{1|0} F'(F P_{1|0} F')^{-1}FP_{1|0}] G' + W $$
