Encoding Drift in State Space Models for ETF Price and NAV
Article Quant Q&A · Author: Grillo
Summary
The document presents a state space formulation for modeling an ETF’s observed log price and net asset value against an unobserved log true value. The true value follows a random walk with drift, while price and NAV depend on their own lags, the true value, and separate observation errors. The author asks how to represent nonzero error means and how to interpret a deterministic constant included in the state vector when estimating the model with R’s dlm package.
The document provides the proposed transition and measurement matrices, an R likelihood function, reported maximum likelihood estimates, and filtered state output. It then works through a simpler random walk with drift to argue that a constant state component should remain one under filtering. The reported fifth filtered component varies, prompting a question about the implementation. These materials illustrate how state augmentation can encode a constant drift, but the document does not resolve the filtering discrepancy or establish whether the estimates are correctly specified. The ETF application is attributed to a published study, though no empirical evaluation is supplied here.
Key ideas
- A random walk with drift can be represented by augmenting the state with a constant component.
- The proposed model relates observed ETF price and NAV to an unobserved true value and their own lagged values.
- The author estimates model parameters with maximum likelihood using the dlm package.
- The filtered constant state component varies in the reported output, despite having deterministic transition and zero process variance.
- The document raises modeling and implementation questions but does not provide definitive answers.
Tags
Full text
# DLM package in R to estimate a state space model with drifts
# 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.
- 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 elsewhere, 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 $$Shown in full with attribution under the source's licence. Licence: CC BY-SA 4.0 (Stack Exchange)
This summary was written by Stratmill's research agent from the original; it is not a copy of the source.