0
$\begingroup$

Here is the problem :-

We have an AR(1) process, $x[t]$,

ie,

$(x[t] - \mu) = \phi(x[t-1]-\mu) + \epsilon_x[t] $

where $Var(\epsilon_x[t]) = \sigma_x^2$ and $Mean(\epsilon_x[t])=0$

ie. $x[t] = (\mu - \phi * \mu) + \phi x[t-1] + \epsilon_x[t] $

ie. $x[t] = \mu(1 - \phi) + \phi x[t-1] + \epsilon_x[t] $

For convenience I call $\mu(1-\phi)$ as $c$.

ie. $x[t] = c + \phi x[t-1] + \epsilon_x[t] $

We observe the 4 period moving average $y[t]$ of $x[t]$ with noise,

ie.

$ y[t] = 1/4*( x[t] + x[t-1] + x[t-2] + x[t-3] ) + \epsilon_y[t]$

where $Var(\epsilon_y[t]) = \sigma_y^2$ and $Mean(\epsilon_y[t])=0$

The problem is to go from $y[t]$ to $x[t]$.

Here is my attempt.

I create the state space formulation of this problem.

The measurement equation being :-

$y[t] = .25 x[t] + .25 x[t-1] + .25 x[t-2] + .25 x[t-3] + \epsilon_y[t]$

The state vector here is $ \begin{bmatrix} x[t] \\ x[t-1] \\ x[t-2] \\ x[t-3] \\ 1 \\ \end{bmatrix} $

The state equation is :-

$ \begin{bmatrix} x[t] \\ x[t-1] \\ x[t-2] \\ x[t-3] \\ 1 \end{bmatrix} = \begin{bmatrix} \phi && 0 && 0 && 0 && c = \mu * (1-\phi) \\ 0 && \phi && 0 && 0 && c = \mu * (1-\phi) \\ 0 && 0 && \phi && 0 && c = \mu * (1-\phi) \\ 0 && 0 && 0 && \phi && c = \mu * (1-\phi) \\ 0 && 0 && 0 && 0 && 1 \\ \end{bmatrix} \begin{bmatrix} x[t-1] \\ x[t-2] \\ x[t-3] \\ x[t-4] \\ 1 \\ \end{bmatrix} + \begin{bmatrix} \epsilon_x[t] \\ \epsilon_x[t-1] \\ \epsilon_x[t-2] \\ \epsilon_x[t-3] \\ 0 \end{bmatrix} $

The variance-covariance matrix of the error matrix above is :-

$ \begin{bmatrix} \sigma_x^2 && 0 && 0 && 0 && 0 \\ 0 && \sigma_x^2 && 0 && 0 && 0 \\ 0 && 0 && \sigma_x^2 && 0 && 0 \\ 0 && 0 && 0 && \sigma_x^2 && 0 \\ 0 && 0 && 0 && 0 && 0 \end{bmatrix} $

In the language of dlm, I have computed GG and W.

$ GG = \begin{bmatrix} \phi && 0 && 0 && 0 && c = \mu * (1-\phi) \\ 0 && \phi && 0 && 0 && c = \mu * (1-\phi) \\ 0 && 0 && \phi && 0 && c = \mu * (1-\phi) \\ 0 && 0 && 0 && \phi && c = \mu * (1-\phi) \\ 0 && 0 && 0 && 0 && 1 \end{bmatrix} $

$ W = \begin{bmatrix} \sigma_x^2 && 0 && 0 && 0 && 0 \\ 0 && \sigma_x^2 && 0 && 0 && 0 \\ 0 && 0 && \sigma_x^2 && 0 && 0 \\ 0 && 0 && 0 && \sigma_x^2 && 0 \\ 0 && 0 && 0 && 0 && 0 \end{bmatrix} $

Next I write the measurement equation in terms of the state vector at time t.

$ y[t] = 1/4*( x[t] + x[t-1] + x[t-2] + x[t-3] ) + \epsilon_y[t]$

ie.

$ y[t] = 1/4 * x[t] + 1/4 * x[t-1] + 1/4 * x[t-2] + 1/4 * x[t-3] + 0 * 1 + \epsilon_y[t]$

ie.

$y[t] = \begin{bmatrix} .25 && .25 && .25 && .25 && 0 \end{bmatrix} * \begin{bmatrix} x[t] \\ x[t-1] \\ x[t-2] \\ x[t-3] \\ 1 \\ \end{bmatrix} + \epsilon_y[t] $

In terms of dlm,

the $FF$ matrix is $ = \begin{bmatrix} .25 && .25 && .25 && .25 && 0 \end{bmatrix} $

and $V$ is $\begin{bmatrix} \sigma_y^2 \end{bmatrix} $

Here is the R program which implements the above

library(dlm)
library(zoo)

# Simulating the data.

set.seed(123)
x <- as.zoo( 6 + arima.sim(model=list(ar=c(.3)), n=1000, sd=2))
y <- zoo::rollapply(x, width=4, FUN=mean, align="right") + rnorm(1000-3, sd=1)
    
# Set parameter restrictions
# parm[1] = Phi, the AR1 parameter of the x series.I have constrained this to be between -1 and 1.
# parm[2] = the error variance in the x series. I have constrained this to be positive.
# parm[3] = the error variance in the y series. I have constrained this to be positive.
# parm[4] = mu is the mean of the x series. I have constrained this to be positive. I know from the physical interpretation of the problem that
# this is positive. In theory this could be negative. It's not a big deal.
# The intercept of the x series is mu*(1-Phi).

parm_rest <- function(parm){
    return(
        c(1-2*exp(parm[1])/(1+exp(parm[1])),
          exp(parm[2]),
          exp(parm[3]),
          exp(parm[4]))
    )
}

# We setup the state space model.

ssm1 <- function(parm){
    
    parm <- parm_rest(parm)

    GG1  = diag(rep(parm[1],4))
    GG1  = cbind(GG1,rep((1-parm[1])*parm[4],4))
    GG1  = rbind(GG1,c(rep(0,4),1))
    W1 = diag(c(rep(parm[2],4),0))
    return(
        dlm(
            FF = matrix(c(rep(.25,4),0),nr=1),
            V = parm[3],
            GG = GG1,
            W = W1,
            m0 = matrix(c(rep(0,4),1),nr=5),
            # I set the first 4 state variables = x[t]=...=x[t-3]=0
            # the 5th state variable = 1
            # C0 = diag(c(rep(solve(1-parm[1]^2)*parm[2],4),0))
            # parm[1] is between -1 and 1 and can be zero. That is why I did NOT do the above as it may lead to division by 0.
            C0 = diag(c(rep(1000,4),0))
            # There is some uncertainty in the x[0],...,x[3].
            # There is NO uncertainty in the last state variable as it is always = 1 
        )
    )
}

# estimate parameters
fit1 <- dlmMLE(y,parm=c(1,1,1,1),build=ssm1,hessian=T)

# get estimates

coef <- parm_rest(fit1$par)

# get standard errors using delta method
dg1 <- -2 * exp(fit1$par[1])/(1+exp(fit1$par[1]))^2
dg2 <- exp(fit1$par[2])
dg3 <- exp(fit1$par[3])
dg4 <- exp(fit1$par[4])
dg <- diag(c(dg1,dg2,dg3,dg4))
var <- dg%*%solve(fit1$hessian)%*%dg

# print results

coef; sqrt(diag(var))

coef + (2 *  sqrt(diag(var)))
coef - (2 *  sqrt(diag(var)))

Here is the output

> coef; sqrt(diag(var))
[1] 0.6532458 4.4895771 0.6795316 6.0745278
[1] 0.04242776 0.66841183 0.12774506 0.10030250
> 
> coef + (2 *  sqrt(diag(var)))
[1] 0.7381013 5.8264007 0.9350217 6.2751328
> coef - (2 *  sqrt(diag(var)))
[1] 0.5683903 3.1527534 0.4240415 5.8739228
> 

The first parameter .3 is NOT between [.56,.73]

The error variance in x = 4 is between [3.15,5.82 ]

The error variance in y = 1 is not between [.42,.93 ]

The mean of x = 6 IS between [5.8,6.2]

This behavior is repeated on repeated sampling. In particular the $\phi$ is not covered in the 2 standard deviations of the coef[1].

Can some tell me where I am mistaken ? Theory or programming? I could have made a mistake but I do not see where I am wrong.

$\endgroup$
0
1
$\begingroup$

The state equation is incorrect. You say: $$ \begin{bmatrix} x[t] \\ x[t-1] \\ x[t-2] \\ x[t-3] \\ 1 \end{bmatrix} = \begin{bmatrix} \phi && 0 && 0 && 0 && c = \mu * (1-\phi) \\ 0 && \phi && 0 && 0 && c = \mu * (1-\phi) \\ 0 && 0 && \phi && 0 && c = \mu * (1-\phi) \\ 0 && 0 && 0 && \phi && c = \mu * (1-\phi) \\ 0 && 0 && 0 && 0 && 1 \\ \end{bmatrix} \begin{bmatrix} x[t-1] \\ x[t-2] \\ x[t-3] \\ x[t-4] \\ 1 \\ \end{bmatrix} + \begin{bmatrix} \epsilon_x[t] \\ \epsilon_x[t-1] \\ \epsilon_x[t-2] \\ \epsilon_x[t-3] \\ 0 \end{bmatrix} $$

If we call the last term term $\varepsilon_t$, your model doesn't know that the first entry of $\varepsilon_t$ should be exactly equal to the 2nd entry of $\varepsilon_{t+1}$, for example. It assumes the $\varepsilon_t$ are independent across time and so will not respect all those implicit dependencies.

The usual way to cast an $\text{AR}(p)$ model in state space form is the "companion form", which in this case would look like this:

$$ \begin{bmatrix} x[t] \\ x[t-1] \\ x[t-2] \\ x[t-3] \\ 1 \end{bmatrix} = \begin{bmatrix} \phi && 0 && 0 && 0 && c \\ 1 && 0 && 0 && 0 && 0 \\ 0 && 1 && 0 && 0 && 0 \\ 0 && 0 && 1 && 0 && 0 \\ 0 && 0 && 0 && 0 && 1 \\ \end{bmatrix} \begin{bmatrix} x[t-1] \\ x[t-2] \\ x[t-3] \\ x[t-4] \\ 1 \\ \end{bmatrix} + \begin{bmatrix} \epsilon_x[t] \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} $$

The first row is then your $\text{AR}(1)$ model and the other rows do the work of "carrying" forward previous values that you need to appear in the observation equation.

As pointed out by @Yves, a better way to handle the constant is then:

$$ \begin{bmatrix} x[t] \\ x[t-1] \\ x[t-2] \\ x[t-3] \\ \mu[t] \end{bmatrix} = \begin{bmatrix} \phi && 0 && 0 && 0 && 1 \\ 1 && 0 && 0 && 0 && 0 \\ 0 && 1 && 0 && 0 && 0 \\ 0 && 0 && 1 && 0 && 0 \\ 0 && 0 && 0 && 0 && 1 \\ \end{bmatrix} \begin{bmatrix} x[t-1] \\ x[t-2] \\ x[t-3] \\ x[t-4] \\ \mu[t-1] \\ \end{bmatrix} + \begin{bmatrix} \epsilon_x[t] \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} $$

$\endgroup$
5
  • 1
    $\begingroup$ Good. I would use $\mu$ as the fourth state component with the transition $\mu[t]= \mu[t-1]$ and a fourth element set to 1 in the observation matrix. The initial variance for $\mu[t]$ should be large or diffuse. Then $\mu$ will no longer be an hyper-parameter (as are $\phi$ and the variances) and the estimate of $\mu$ would simply be the filtered $\widehat{\mu}[T|T]$ where $T$ is the last time, provided that the ML estimation of the hyperparameters have converged. $\endgroup$ – Yves May 19 at 12:15
  • $\begingroup$ Assuming the the observation remains with a zero fourth element I think that the fourth component of the answer is $c[t]$ rather than $\mu[t]$. An equivalent change is to drop the one in the first line of the transition and use a one in the observation matrix (4th position). The first 3 state elements should then fluctuate around zero. $\endgroup$ – Yves May 19 at 12:26
  • $\begingroup$ Dear Yves and Chris, do you mean to use $ c = \mu ( 1- \phi ) $ as the coefficient corresponding to the fifth component which is 1 ? If I do that the standard error in c and phi will be known ( it will be spit out by the program ) , but how will you compute the standard error in $ \mu $ ? That is why I explicitly modelled $ \mu, \phi $ in place of $c$. What do you think about this ? $\endgroup$ – user2338823 May 19 at 12:32
  • $\begingroup$ Sorry typo in my comments: "fourth" is to be read as "fifth". If $\mu$ is known to be zero then the state would be 4-dimensional with a transition in "companion form" as in @Chris Haug answer. If a constant $\mu$ is to be added to $x[t]$ hence to $y[t]$ we can make a "state augmentation" using $\mu$ as a new state element. The filltered $\widehat{\mu}(t | t)$ will converge to the estimate of $\mu$ (for given hyper parameters). The standard error can be derived from the filtered variance at $T$. You could use $c$ in the state to shift $x[t]$ as in the answer. This is profile-likelihood. $\endgroup$ – Yves May 19 at 13:02
  • $\begingroup$ Hi again, I have a small query. We could model $\mu$ as hyper parameter as in the original query OR we could make $\mu[t]$ the last element of the state vector as suggested by Yves. My query is : why is modelling $\mu[t]$ as a state variable a better idea? In both cases I would get the mean and the standard error of $\mu$. Why is the latter approach preferable? $\endgroup$ – user2338823 May 22 at 3:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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