0
$\begingroup$

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:

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 $$

$\endgroup$
15
  • $\begingroup$ I don't have time to look at the R code carefully but, if you're error terms have non-zero means then you need to specify the means somewhere in the formulation ( I think the FF matrix would be a good location ) . I don't see that being done somewhere but maybe you mean something else when you say "independent but with a drift". Do you have Giovanni's dlm UseR book ? It's chock full of good examples. $\endgroup$
    – mlofton
    Commented Nov 12, 2020 at 20:05
  • $\begingroup$ Thank you for your suggestion: I have already done that: The two means of the measurement equation variables are indeed in FF and the mean of the state variable is in GG. These means are added to the measurement variables and the state variable by multiplying the matrices with the state vector entry in [5,1] that is equal to 1. Of course, this value should remain to be equal to 1 but somehow it changes based on the output from dlmFilter(). I don’t understand why, which led me to ask my question #2. Maybe I don’t understand dlmFilter() correctly. $\endgroup$
    – Grillo
    Commented Nov 12, 2020 at 21:15
  • $\begingroup$ Hi: I'll print it out and try to read during my travels today. It's been a long time but I used to do this sort of thing a lot. Can you get your hands on giovanni's dlm useR book ? It's very nice book and that may help a lot. I'll be in touch but I can't guarantee that I can help you. $\endgroup$
    – mlofton
    Commented Nov 13, 2020 at 12:06
  • $\begingroup$ Can you write out the observation equation and system equation without matrices because they don't feel correct. F is like the design matrix in regression ( the betas ), $\theta$ is the state vector., $G$ is a matrix representing how the state vector evolves. So, you need to define $F$ and $G$ and $\theta$ correctly. Then, once they are defined, things become easier. Note that the error terms are not state variables but you have $v_t$ evolving as if it's a state variable. ( with $v_{t-1}$ fixed ). Is $v_t$ a state variable or an error term ? Also, where is the model coming from ? $\endgroup$
    – mlofton
    Commented Nov 14, 2020 at 1:26
  • $\begingroup$ my DLM useR book is not where I currently am but i found a pdf of it at the link below. pages 41-46 of it have some examples of some basic models that still might be helpful. I would take a look at those and also write your model out without matrices first. Then, if you can send the model without matrices to me, that would be good. researchgate.net/publication/… $\endgroup$
    – mlofton
    Commented Nov 14, 2020 at 1:40

1 Answer 1

0
$\begingroup$
  1. As for constant means, I think you would need a fixed term in the state vector.

  2. When you set variance zero for an element of the state vector, you make it fixed; but its estimate will evolve over time, as more observations are used. Think of a simple regression with only an intercept (fixed parameter): the estimate will change as you include more observations.

What is true, I think, is that the smoothed estimate would be constant (as it depends on the whole sample at any point in time).

$\endgroup$
13
  • $\begingroup$ Thank you for your answer. I appreciated your help but I do not believe this is correct. The filtered value of the constant "1" in the last position of the state vector should remain the same value "1" in the filtered state vector $v_{t|t}$ for all $t$. I am happy to write this out. Alternatively, I can post the filtered values based on the MLE estimates above and the model in dlm(). Please let me know. Alternatively, I can show you the manually updated inference about state vector, which I beleive to be the output of dlmFilter(). $\endgroup$
    – Grillo
    Commented Nov 17, 2020 at 6:12
  • $\begingroup$ I am afraid I don't understand you, then. You already posted the filtered values, and found X5 varying, which apparently prompted your question. Now you say this is not the case? $\endgroup$
    – F. Tusell
    Commented Nov 17, 2020 at 9:59
  • $\begingroup$ The values I posted more recently are the manually filtered values of the state vector using the updating procedure in Hamilton. They are different than what dlmFilter() produces. In particular they retain the constant value of "1" in the last entry of the filtered state vector. I have added some clarification using a simpler model to show (for the first few steps) that using the recursive process to obtain filtered values of the state vector that that the entry of a "1" should be preserved in the filtered state vector. $\endgroup$
    – Grillo
    Commented Nov 17, 2020 at 19:05
  • $\begingroup$ In general (and this is the assumption of dlm unless you override it) we start the Kalman filter recursion with a diffuse prior. In dlm this is approximated by an initial state covariance matrix diagonal with terms "large" $(10^7)$, which allow for the state vector elements to change even if the corresponding noises have variance zero. $\endgroup$
    – F. Tusell
    Commented Nov 18, 2020 at 10:03
  • $\begingroup$ Thank you for your reply. I obviously know that the "1" entry in the state vector needs to remain "1" if the model, as set up, were to make sense. Since you are suggesting that this is not feasible, for example by setting the prior for all but that one element to 10E7 and to zero for that element (which I have done with no difference in the result), do you have a suggestion how to specify the state-space model allowing for a constant mean in the latent variable ($v_t$) in my case? $\endgroup$
    – Grillo
    Commented Nov 18, 2020 at 16:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.