Filtering data with the Kalman Filter

Goals

Today we will look at the GAUSS kalman filter procedure, which is included in the TSMT application module.

In this tutorial we will use an AR(2) example to examine

1. State space representation
2. How to define system matrices
3. How to use the GAUSS kalmanFilter procedure
4. How to interpret the results of the kalmanResult structure.

Introduction

The Kalman filter is a tool that allows us to determine the optimal estimates of an unobserved state vector, $\alpha_t$, using the observed information at time $t$.

Today will demonstrate the use of the Kalman filter to predict an AR(2) model from simulated data. We generate our sample AR(2) model using the data generating process

$$y_t = 0.3 y_{t-1} - 0.6 y_{t-2} + \epsilon_t .$$

State space representation

The GAUSS kalmanFilter procedure requires that models take a state space representation. Therefore, we must understand how to translate our problem into proper state space representation.

State space representation is composed of three pieces:

• An observation equation.
• A measurement equation.
• An initial state distributions.

The observation equation represents an observed time-series vector

$$y_t = d_t + Z_t \alpha_t + \epsilon_t\\ \epsilon_t \sim N(0, H_t)$$

The transition equation represents an unobserved state vector

$$\alpha_{t+1} = c_t + T_t \alpha_t + R\eta_t\\ \eta_t \sim N(0, Q_t)$$

Finally, the initial state distribution represents starting values of the state vector. It includes the initial prior state mean, $a_0$, and the initial prior state covariance, $P_0$.

$$\alpha_0 \sim N(a_0, P_0)$$

State space representation of an AR(2) model

Putting models into state space representation is somewhat of an art form and there can be many different representations of the same model.

Let's consider putting a general AR(2) model into the state space representation starting from the standard representation

$$y_t = \phi_1 y_{t-1} + \phi_2 y_{t-2} + \epsilon_t .$$

AR(2) representation #1

For our first representation we will define

$$\alpha_t = \begin{bmatrix} y_t\\ y_{t-1}\end{bmatrix}$$

then our measurement equation becomes

$$y_t = \begin{bmatrix} 1 & 0\end{bmatrix} \alpha_t$$

and our transition equation becomes

$$ \alpha_t = \begin{bmatrix} \phi_1 & \phi_2 \\ 1 & 0 \\ \end{bmatrix} \alpha_{t-1} + \begin{bmatrix} 1\\ 0\end{bmatrix} \eta_t $$

To see how this works, let's look more closely at the matrix multiplication. Consider first plugging $\alpha_t$ into our measurement equation

$$y_t = \begin{bmatrix} 1 & 0\end{bmatrix} \begin{bmatrix} y_t\\ y_{t-1}\end{bmatrix} = 1*y_t + 0 * y_{t-1} = y_t$$

and

$$ \alpha_t = \begin{bmatrix} \phi_1 & \phi_2 \\ 1 & 0 \\ \end{bmatrix} \alpha_{t-1} + \begin{bmatrix} 1\\ 0\end{bmatrix} \eta_t =\\ \begin{bmatrix} \phi_1 & \phi_2 \\ 1 & 0 \\ \end{bmatrix} \begin{bmatrix} y_{t-1}\\ y_{t-2}\end{bmatrix} + \begin{bmatrix} 1\\ 0\end{bmatrix} \eta_t =\\ \begin{bmatrix} \phi_1 y_{t-1} + \phi_2 y_{t-2}\\ 1*y_{t-1} + 0*y_{t-2}\end{bmatrix} + \begin{bmatrix} 1\\ 0\end{bmatrix} \eta_t =\\ \begin{bmatrix} y_t\\ y_{t-1}\end{bmatrix} $$

which simplifies to

$$y_t = \phi_1 y_{t-1} + \phi_2 y_{t-2} + \eta_t\\ y_{t-1} = y_{t-1}$$

AR(2) representation #2

We could also an alternative representation where we let

$$\alpha_t = \begin{bmatrix} y_t\\ \phi_2 y_{t-1}\end{bmatrix}$$

then our measurement equation becomes

$$y_t = \begin{bmatrix} 1 & 0\end{bmatrix} \alpha_t$$

and our transition equation becomes

$$ \alpha_t = \begin{bmatrix} \phi_1 & 1 \\ \phi_2 & 0 \\ \end{bmatrix} \alpha_{t-1} + \begin{bmatrix} 1\\ 0\end{bmatrix} \eta_t $$

Using the GAUSS Kalman filter

The GAUSS kalmanFilter procedure takes 10 inputs. The first of these inputs is the observed data, $y_t$. We will use the TSMT application module procedure simarmamt to generate our data.

Observed AR(2) data

library tsmt;

// Simulate AR(2) data
seed = 527453;

// Intercept
b = 0.3|-0.6;

// AR terms
p = 2;

// MA terms
q = 0;

// Constant
const = 0;

// Trend
trend = 0;

// Number of observations
n = 150;

// Number of replications
k = 1;

// Standard deviation of error process
std = 1;

// Generate data
y = simarmamt(b, p, q, const, trend, n, k, std, seed)';

System matrices

The 2nd through 8th inputs into the kalmanFilter procedure are based on the system matrices of the general state space representation.

GAUSS kalmanFilter system matrices inputs
$Z$<\td>	The design matrix.
$d$<\td>	The measurement equation intercept.
$H$<\td>	The measurement equation disturbance covariance matrix.
$T$<\td>	The transition matrix.
$c$<\td>	The state transition equation intercept.
$R$<\td>	The selection matrix.
$Q$<\td>	The state disturbance covariance matrix.

/*
** Use the two true
** phi values
*/
phi = 0.3|-0.6;

// State disturbance sigma
sigma2_q = 1;

// Observation disturbance sigma
sigma2_h = 1;

// Design matrix
Z = {1 0};

// Measurement equation intercept
d = 0;

// Measurement equation disturbance
H = sigma2_h;

// Transition matrix
T = {phi[1] phi[2], 1 0};

// Transition equation intercept
c = 0;

// Transition equation disturbance
Q = sigma2_q;

// Selection matrix
R = {1, 0};

Initial state distribution

Before we can run the Kalman filter we must initialize the state vector. This requires specifying both $a_0$, the initial prior state mean, and $P_0$, the initial prior state covariance.

// Initial state mean
a_0 = 0|0;

// Initial state covariance
P_0 = {1 0, 0 1};

Calling the kalmanFilter procedure

The kalmanFilter returns output to a kalmanResult structure whose members will be discussed in the next section. The important thing to note now is that before calling the procedure we must declare an instance of the kalmanResult structure.

struct kalmanResult rslt;
rslt = kalmanFilter(yt, Z, d, H, T, c, R, Q, a_0, p_0);

Summary of how to use kalmanFilter

We have covered a lot of information in this section. It may be helpful to have a summary of the steps we have to take to use the kalmanFilter procedure:

1. Determine state space representation.
2. Define the observed data vector.
3. Define the system matrices.
4. Define the initial state mean and initial state covariance.
5. Declare the kalmanResult structure and call the kalmanFilter procedure.

The kalmanResult structure

All output from the kalmanFilter procedure is stored in the kalmanResult structure. The members of the kalmanResult instance kout include:

kout.filtered_state	The filtered states.
kout.filtered_state_cov	Filtered state covariances.
kout.predicted_state	The predicted states.
kout.predicted_state_cov	Predicted state covariances.
kout.forecast	Forecasts.
kout.forecast_error	Forecast error.
kout.forecast_error_cov	Forecast error covariances.
kout.loglikelihood	Computed log likelihood.

Let's plot the results stored in kout.predicted_state against a scatter plot of the original data (view the code to produce the graphs here).

We can see that the predicted state closely follows our observed data.

Conclusion

Congratulations! You have completed the Kalman filter tutorial. In this tutorial, we used an AR(2) example to examine

1. State space representation
2. How to define system matrices
3. How to use the GAUSS kalmanFilter procedure
4. The members of the kalmanResult structure.

For your convenience, the entire code is available below.

library tsmt;

// Simulate AR(2) data
seed = 527453;

// Intercept
b = 0.3|-0.6;

// AR terms
p = 2;

// MA terms
q = 0;

// Constant
const = 0;

// Trend
trend = 0;

// Number of observations
n = 150;

// Number of replications
k = 1;

// Standard deviation of error process
std = 1;

// Generate data
y = simarmamt(b, p, q, const, trend, n, k, std, seed)';

/*
** Use the two true
** phi values
*/
phi = 0.3|-0.6;

// State disturbance sigma
sigma2_q = 1;

// Observation disturbance sigma
sigma2_h = 1;

// Design matrix
Z = {1 0};

// Measurement equation intercept
d = 0;

// Measurement equation disturbance
H = sigma2_h;

// Transition matrix
T = {phi[1] phi[2], 1 0};

// Transition equation intercept
c = 0;

// Transition equation disturbance
Q = sigma2_q;

// Selection matrix
R = {1, 0};

// Initial state mean
a_0 = 0|0;

// Initial state covariance
P_0 = {1 0, 0 1};

struct kalmanResult rslt;
rslt = kalmanFilter(yt, Z, d, H, T, c, R, Q, a_0, p_0);