# Using a state space model to invert a moving average

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