Chapter 4 Stan Functions

State space functionality for Stan is provided as a set of user-defined functions.

Add the following line to the Stan model file in which depends on these functions.

functions {
  #include ssm.stan
  // other functions ...
}

To actually include the functions in the model, you need to use the function stanc_builder, instead of stan or stanc:

model <- stanc_builder("yourmodel.stan", isystem = "path/to/ssm/")
stan(model_code = model$model_code)

4.1 Utility Functions

4.1.1 to_symmetric_matrix

Parameters:

• x An \(n \times n\) matrix.

Return Value: An \(n \times n\) symmetric matrix: \(0.5 (x + x')\).

Ensure a matrix is symmetrix.

matrix to_symmetric_matrix(matrix x) {
  return 0.5 * (x + x ');
}

4.1.2 to_matrix_colwise

Parameters:

• vector v An \(n \times m\) vector.
• int m Number of rows in the vector
• int n Number of columns in the vector

Return Value: matrix A \(m \times n\) matrix containting the elements from v

Convert vector to a matrix (column-major).

matrix to_matrix_colwise(vector v, int m, int n) {
  matrix[m, n] res;
  for (j in 1:n) {
    for (i in 1:m) {
      res[i, j] = v[(j - 1) * m + m];
    }
  }
  return res;
}

4.1.3 to_matrix_rowwise

Parameters:

• vector v An \(n \times m\) vector.
• int m Number of rows in the matrix.
• int n Number of columns in the matrix.

Return Value: matrix A \(m \times n\) matrix containting the elements from v

Convert vector to a matrix (row-major).

matrix to_matrix_rowwise(vector v, int m, int n) {
  matrix[m, n] res;
  for (i in 1:n) {
    for (j in 1:m) {
      res[i, j] = v[(i - 1) * n + n];
    }
  }
  return res;
}

4.1.4 to_vector_colwise

Parameters:

• matrix x An \(n \times m\) matrix.

Return Value: vector with \(n m\) elements.

Convert a matrix to a vector (column-major)

vector to_vector_colwise(matrix x) {
  vector[num_elements(x)] res;
  int n;
  int m;
  n = rows(x);
  m = cols(x);
  for (i in 1:n) {
    for (j in 1:m) {
      res[n * (j - 1) + i] = x[i, j];
    }
  }
  return res;
}

4.1.5 to_vector_rowwise

Parameters:

• matrix x An \(n \times m\) matrix.

Return Value: vector with \(n m\) elements.

Convert a matrix to a vector (row-major)

vector to_vector_rowwise(matrix x) {
  vector[num_elements(x)] res;
  int n;
  int m;
  n = rows(x);
  m = cols(x);
  for (i in 1:rows(x)) {
    for (j in 1:cols(x)) {
      res[(i - 1) * m + j] = x[i, j];
    }
  }
  return res;
}

4.1.6 symmat_size

Parameters:

• matrix x An \(m \times m\) matrix.

Return Value: int The number of unique elements

Calculate the number of unique elements in a symmetric matrix

The number of unique elements in an \(m \times m\) matrix is \((m \times (m + 1)) / 2\).

int symmat_size(int n) {
  int sz;
  // This calculates it iteratively because Stan gives a warning
  // with integer division.
  sz = 0;
  for (i in 1:n) {
    sz = sz + i;
  }
  return sz;
}

4.1.7 find_symmat_dim

Parameters:

• int n The number of unique elements in a symmetric matrix.

Return Value: int The dimension of the associated symmetric matrix.

Given vector with \(n\) elements containing the \(m (m + 1) / 2\) elements of a symmetric matrix, return \(m\).

int find_symmat_dim(int n) {
  // This could be solved by finding the positive root of $m = m (m + 1)/2 but
  // Stan doesn't support all the functions necessary to do this.
  int i;
  int remainder;
  i = 0;
  while (n > 0) {
    i = i + 1;
    remainder = remainder - i;
  }
  return i;
}

4.1.8 vector_to_symmat

Parameters:

• vector x The vector with the unique elements
• int n The dimensions of the returned matrix: \(n \times n\).

Return Value: matrix An \(n \times n\) symmetric matrix.

Convert a vector to a symmetric matrix

matrix vector_to_symmat(vector x, int n) {
  matrix[n, n] m;
  int k;
  k = 1;
  for (j in 1:n) {
    for (i in 1:j) {
      m[i, j] = x[k];
      if (i != j) {
        m[j, i] = m[i, j];
      }
      k = k + 1;
    }
  }
  return m;
}

4.1.9 symmat_to_vector

Parameters:

• vector x An \(n \times n\) matrix.

Return Value: vector A \(n (n + 1) / 2\) vector with the unique elements in \(x\).

Convert an \(n \times n\) symmetric matrix to a length \(n (n + 1) / 2\) vector containing its unique elements.

vector symmat_to_vector(matrix x) {
  vector[symmat_size(rows(x))] v;
  int k;
  k = 1;
  // if x is m x n symmetric, then this will return
  // only parts of an m x m matrix.
  for (j in 1:rows(x)) {
    for (i in 1:j) {
      v[k] = x[i, j];
      k = k + 1;
    }
  }
  return v;
}

4.2 Filtering

Functions used in filtering and log-likelihood calculations.

4.2.1 ssm_filter_update_a

Parameters:

• vector a An \(m \times 1\) vector with the prected state, \(\vec{a}_t\).
• vector c An \(m \times 1\) vector with the system intercept, \(\vec{c}_t\)
• matrix T An \(m \times m\) matrix with the transition matrix, \(\mat{T}_t\).
• vector v A \(p \times 1\) vector with the forecast error, \(\vec{v}_t\).
• matrix K An \(m \times p\) matrix with the Kalman gain, \(\mat{K}_t\).

Return Value: vector A \(m \times 1\) vector with the predicted state at \(t + 1\), \(\vec{a}_{t + 1}\).

Update the expected value of the predicted state, \(\vec{a}_{t + 1} = \E(\vec{\alpha}_{t + 1} | \vec{y}_{1:t})\),

The predicted state \(\vec{a}_{t + 1}\) is, \[ \vec{a}_{t + 1} = \mat{T}_t \vec{a}_t + \mat{K}_t \vec{v}_t + \vec{c}_t . \]

vector ssm_filter_update_a(vector a, vector c, matrix T, vector v, matrix K) {
  vector[num_elements(a)] a_new;
  a_new = T * a + K * v + c;
  return a_new;
}

4.2.2 ssm_filter_update_P

Parameters:

• matrix P An \(m \times m\) vector with the variance of the prected state, \(\mat{P}_t\).
• matrix Z A \(p \times m\) matrix with the design matrix, \(\mat{Z}_t\).
• matrix T An \(m \times m\) matrix with the transition matrix, \(\mat{T}_t\).
• matrix RQR A \(m \times m\) matrix with the system covariance matrix, \(\mat{R}_t \mat{Q}_t \mat{R}_t'\).
• matrix K An \(m \times p\) matrix with the Kalman gain, \(\mat{K}_t\).

Return Value: matrix An \(m \times 1\) vector with the predicted state at \(t + 1\), \(\vec{a}_{t + 1}\).

Update the expected value of the predicted state, \(\mat{P}_{t + 1} = \Var(\alpha_{t + 1} | \vec{y}_{1:t})\),

The predicted state variance \(\mat{P}_{t + 1}\) is, \[ \mat{P}_{t + 1} = \mat{T}_t \mat{P}_t (\mat{T}_t - \mat{K}_t \mat{Z}_t)' + \mat{R}_t \mat{Q}_t \mat{R}_t' . \]

matrix ssm_filter_update_P(matrix P, matrix Z, matrix T,
                           matrix RQR, matrix K) {
  matrix[rows(P), cols(P)] P_new;
  P_new = to_symmetric_matrix(T * P * (T - K * Z)' + RQR);
  return P_new;
}

4.2.3 ssm_filter_update_v

Parameters:

• matrix P An \(m \times m\) vector with the variance of the prected state, \(\mat{P}_t\).
• matrix Z A \(p \times m\) matrix with the design matrix, \(\mat{Z}_t\).
• matrix T An \(m \times m\) matrix with the transition matrix, \(\mat{T}_t\).
• matrix RQR An \(m \times m\) matrix with the system covariance matrix, \(\mat{R}_t \mat{Q}_t \mat{R}_t'\).
• matrix K An \(m \times p\) matrix with the Kalman gain, \(\mat{K}_t\).

Return Value: vector An \(m \times 1\) vector with the predicted state at \(t + 1\), \(\vec{a}_{t + 1}\).

Update the forcast error, \(\vec{v}_t = \vec{y}_t - \E(\vec{y}_t | \vec{y_{1:(t - 1)}})\)

The forecast error \(\vec{v}_t\) is \[ \vec{v}_t =\vec{y}_t - \mat{Z}_t \vec{a}_t - \vec{d}_t . \]

vector ssm_filter_update_v(vector y, vector a, vector d, matrix Z) {
  vector[num_elements(y)] v;
  v = y - Z * a - d;
  return v;
}

4.2.4 ssm_filter_update_F

Parameters:

• matrix P An \(m \times m\) vector with the variance of the prected state, \(\mat{P}_t\).
• matrix Z A \(p \times m\) matrix with the design matrix, \(\mat{Z}_t\).
• matrix H A \(p \times p\) matrix with the observation covariance matrix, \(\mat{H}_t\).

Return Value: matrix A \(p \times p\) vector with \(\mat{F}_t\).

Update the variance of the forcast error, \(\mat{F}_t = \Var(\vec{y}_t - \E(\vec{y}_t | \vec{y_{1:(t - 1)}}))\)

The variance of the forecast error \(\mat{F}_t\) is \[ \mat{F}_t = \mat{Z}_t \mat{P}_t \mat{Z}_t + \mat{H}_t . \]

matrix ssm_filter_update_F(matrix P, matrix Z, matrix H) {
  matrix[rows(H), cols(H)] F;
  F = quad_form(P, Z') + H;
  return F;
}

4.2.5 ssm_filter_update_Finv

Parameters:

• matrix P An \(m \times m\) vector with the variance of the prected state, \(\mat{P}_t\).
• matrix Z A \(p \times m\) matrix with the design matrix, \(\mat{Z}_t\).
• matrix H A \(p \times p\) matrix with the observation covariance matrix, \(\mat{H}_t\).

Return Value: matrix A \(p \times p\) vector with \(\mat{F}^{-1}_t\).

Update the precision of the forcast error, \(\mat{F}^{-1}_t = \Var(\vec{y}_t - \E(\vec{y}_t | \vec{y_{1:(t - 1)}}))^{-1}\)

This is the inverse of \(\mat{F}_t\).

matrix ssm_filter_update_Finv(matrix P, matrix Z, matrix H) {
  matrix[rows(H), cols(H)] Finv;
  Finv = inverse(ssm_filter_update_F(P, Z, H));
  return Finv;
}

4.2.6 ssm_filter_update_K

Parameters:

• matrix P An \(m \times m\) vector with the variance of the prected state, \(P_t\).
• matrix Z A \(p \times m\) matrix with the design matrix, \(\mat{Z}_t\).
• matrix T An \(m \times m\) matrix with the transition matrix, \(\mat{T}_t\).
• matrix Finv A \(p \times p\) matrix

Return Value: matrix An \(m \times p\) matrix with the Kalman gain, \(\mat{K}_t\).

Update the Kalman gain, \(\mat{K}_t\).

The Kalman gain is \[ \mat{K}_t = \mat{T}_t \mat{P}_t \mat{Z}_t' \mat{F}^{-1}_t . \]

matrix ssm_filter_update_K(matrix P, matrix Z, matrix T, matrix Finv) {
  matrix[cols(Z), rows(Z)] K;
  K = T * P * Z' * Finv;
  return K;
}

4.2.7 ssm_filter_update_L

Parameters:

• matrix Z A \(p \times m\) matrix with the design matrix, \(\mat{Z}_t\)
• matrix T An \(m \times m\) matrix with the transition matrix, \(\mat{T}_t\).
• matrix K An \(m \times p\) matrix with the Kalman gain, \(\mat{K}_t\).

Return Value: matrix An \(m \times m\) matrix, \(\mat{L}_t\).

Update \(L_t\)

\[ \mat{L}_t = \mat{T}_t - \mat{K}_t \mat{Z}_t . \]

matrix ssm_filter_update_L(matrix Z, matrix T, matrix K) {
  matrix[rows(T), cols(T)] L;
  L = T - K * Z;
  return L;
}

4.2.8 ssm_filter_update_ll

Parameters:

• vector v A \(p \times 1\) matrix with the forecast error, \(\vec{v}_t\).
• matrix Finv A \(p \times p\) matrix with variance of the forecast error, \(\mat{F}^{-1}_t\).

Return Value: real An \(m \times m\) matrix, \(L_t\).

Calculate the log-likelihood for a period

The log-likehood of a single observation in a state-space model is \[ \ell_t = - \frac{1}{2} p \log(2 \pi) - \frac{1}{2} \left(\log|\mat{F}_t| + \vec{v}_t' \mat{F}^{-1}_t \vec{v}_t \right) \]

real ssm_filter_update_ll(vector v, matrix Finv) {
  real ll;
  int p;
  p = num_elements(v);
  // det(A^{-1}) = 1 / det(A) -> log det(A^{-1}) = - log det(A)
  ll = (- 0.5 *
        (p * log(2 * pi())
         - log_determinant(Finv)
         + quad_form(Finv, v)
       ));
  return ll;
}

4.3 Filtering

4.3.1 ssm_filter_idx

Parameters:

• int m The number of states
• int p The size of the observation vector \(\vec{y}_t\).

Return Value: int[,] A \(6 \times 3\) integer array containing the indexes of the return values of the Kalman filter.

Indexes of the return values of the Kalman filter functions: ssm_filter.

ssm_filter_idx returns a \(6 \times 3\) integer array with the (length, start index, stop index) of (\(\ell_t\), \(\vec{v}\), \(\vec{F}^-1\), \(\mat{K}\), \(\vec{a}\), \(\mat{P}\)).

value	length	start	stop
\(\ell_t\)	\(1\)	\(1\)	\(1\)
\(\vec{v}\)	\(p\)	\(2\)	\(1 + p\)
\(\mat{F}^{-1}\)	\(p (p + 1) / 2\)	\(2 + p\)	\(1 + p + p (p + 1) / 2\)
\(\mat{K}\)	\(mp\)	\(2 + p + p (p + 1) / 2\)	\(1 + p + p (p + 1) / 2 + mp\)
\(\vec{a}_t\)	\(m\)	\(2 + p + p (p + 1) / 2 + mp\)	\(1 + p + p (p + 1) / 2 + mp + m\)
\(\mat{P}^t\)	\(m (m + 1) / 2\)	\(2 + p + p (p + 1) / 2 + mp + m\)	\(1 + p + p (p + 1) / 2 + mp + m (m + 1) / 2\)

int[,] ssm_filter_idx(int m, int p) {
  int sz[6, 3];
  // loglike
  sz[1, 1] = 1;
  // v
  sz[2, 1] = p;
  // Finv
  sz[3, 1] = symmat_size(p);
  // K
  sz[4, 1] = m * p;
  // a
  sz[5, 1] = m;
  // P
  sz[6, 1] = symmat_size(m);
  // Fill in start and stop points
  sz[1, 2] = 1;
  sz[1, 3] = sz[1, 2] + sz[1, 1] - 1;
  for (i in 2:6) {
    sz[i, 2] = sz[i - 1, 3] + 1;
    sz[i, 3] = sz[i, 2] + sz[i, 1] - 1;
  }
  return sz;
}

4.3.2 ssm_filter_size

Parameters:

• int m The number of states
• int p The size of the observation vector \(\vec{y}_t\).

Return Value: int The number of elements in the vector.

Number of elements in vector containing filter results

int ssm_filter_size(int m, int p) {
  int sz;
  int idx[6, 3];
  idx = ssm_filter_idx(m, p);
  sz = idx[6, 3];
  return sz;
}

4.3.3 ssm_filter_get_loglik

Parameters:

• vector A vector with results from ssm_filter.
• int m The number of states
• int p The size of the observation vector \(\vec{y}_t\).

Return Value: real The log-likelihood \(\ell_t\)

Get the log-likehood from the results of ssm_filter.

real ssm_filter_get_loglik(vector x, int m, int p) {
  real y;
  y = x[1];
  return y;
}

4.3.4 ssm_filter_get_v

Parameters:

• vector A vector with results from ssm_filter.
• int m The number of states
• int p The size of the observation vector \(\vec{y}_t\).

Return Value: vector A \(p \times 1\) vector with the forecast error, \(\vec{v}_t\).

Get the forecast error from the results of ssm_filter.

vector ssm_filter_get_v(vector x, int m, int p) {
  vector[p] y;
  int idx[6, 3];
  idx = ssm_filter_idx(m, p);
  y = segment(x, idx[2, 2], idx[2, 3]);
  return y;
}

4.3.5 ssm_filter_get_Finv

Parameters:

• vector A vector with results from ssm_filter.
• int m The number of states
• int p The size of the observation vector \(\vec{y}_t\).

Return Value: matrix A \(p \times p\) matrix with the forecast precision, \(\mat{F}^{-1}_t\).

Get the forecast precision from the results of ssm_filter.

matrix ssm_filter_get_Finv(vector x, int m, int p) {
  matrix[p, p] y;
  int idx[6, 3];
  idx = ssm_filter_idx(m, p);
  y = vector_to_symmat(segment(x, idx[3, 2], idx[3, 3]), p);
  return y;
}

4.3.6 ssm_filter_get_K

Parameters:

• vector A vector with results from ssm_filter.
• int m The number of states
• int p The size of the observation vector \(\vec{y}_t\).

Return Value: matrix A \(m \times p\) matrix with the Kalman gain, \(\mat{F}^{-1}_t\).

Get the Kalman gain from the results of ssm_filter.

matrix ssm_filter_get_K(vector x, int m, int p) {
  matrix[m, p] y;
  int idx[6, 3];
  idx = ssm_filter_idx(m, p);
  y = to_matrix_colwise(segment(x, idx[4, 2], idx[4, 3]), m, p);
  return y;
}

4.3.7 ssm_filter_get_a

Parameters:

• vector A vector with results from ssm_filter.
• int m The number of states
• int p The size of the observation vector \(\vec{y}_t\).

Return Value: vector An \(m \times 1\) vector with the expected value of the predicted state, \(\E(\vec{\alpha}_t | \vec{y}_{1:(t-1)}) = \vec{a}_t\).

Get the expected value of the predicted state from the results of ssm_filter.

vector ssm_filter_get_a(vector x, int m, int p) {
  vector[m] y;
  int idx[6, 3];
  idx = ssm_filter_idx(m, p);
  y = segment(x, idx[5, 2], idx[5, 3]);
  return y;
}

4.3.8 ssm_filter_get_P

Parameters:

• vector A vector with results from ssm_filter.
• int m The number of states
• int p The size of the observation vector \(\vec{y}_t\).

Return Value: matrix An \(m \times m\) matrix with the variance of the predicted state, \(\Var(\vec{\alpha}_t | \vec{y}_{1:(t-1)}) = \mat{P}_t\).

Get the variance of the predicted state from the results of ssm_filter.

matrix ssm_filter_get_P(vector x, int m, int p) {
  matrix[m, m] y;
  int idx[6, 3];
  idx = ssm_filter_idx(m, p);
  y = vector_to_symmat(segment(x, idx[6, 2], idx[6, 3]), m);
  return y;
}

4.3.9 ssm_filter

Parameters:

• vector[] y Observations, \(\vec{y}_t\). An array of size \(n\) of \(p \times 1\) vectors.
• vector[] d Observation intercept, \(\vec{d}_t\). An array of \(p \times 1\) vectors.
• matrix[] Z Design matrix, \(\mat{Z}_t\). An array of \(p \times m\) matrices.
• matrix[] H Observation covariance matrix, \(\mat{H}_t\). An array of \(p \times p\) matrices.
• vector[] c State intercept, \(\vec{c}_t\). An array of \(m \times 1\) vectors.
• matrix[] T Transition matrix, \(\mat{T}_t\). An array of \(m \times m\) matrices.
• matrix[] R State covariance selection matrix, \(\mat{R} _t\). An array of \(p \times q\) matrices.
• matrix[] Q State covariance matrix, \(\mat{Q}_t\). An array of \(q \times q\) matrices.
• vector a1 Expected value of the intial state, \(a_1 = \E(\alpha_1)\). An \(m \times 1\) matrix.
• matrix P1 Variance of the initial state, \(P_1 = \Var(\alpha_1)\). An \(m \times m\) matrix.

Return Value: vector[] Array of size \(n\) of \((1 + p + p (p + 1) / 2 + mp + m + m (m + 1) / 2) \times 1\) vectors in the format described in ssm_filter_idx.

Kalman filter

For d, Z, H, c, T, R, Q the array can have a size of 1, if it is not time-varying, or a size of \(n\) (for d, Z, H) or \(n - 1\) (for c, T, R, Q) if it is time varying.

ssm_filter runs a forward filter on the state space model and calculates,

• log-likelihood for each observation, \(\ell_t\).
• Forecast error, \(\vec{v}_t = \vec{y}_t - \E(\vec{y}_t | \vec{y}_{1:(t -1)})\).
• Forecast precision, \(\mat{F}^{-1}_t\).
• Kalman gain, \(\mat{K}_t\).
• Predicted states, \(\vec{a}_t = \E(\vec{\alpha}_t | \vec{y}_{1:(t -1)})\).
• Variance of the predicted states, \(\mat{P}_t = \Var(\vec{\alpha}_t | \vec{y}_{1:(t -1)})\).

The results of Kalman filter for a given are returned as a \(1 + p + p (p + 1) / 2 + m p + m (m + 1) / 2\) vector for each time period, where \[ (\ell_t, \vec{v}_t', \VEC(\mat{F}^{-1}_t)', \VEC(\mat{K}_t)', \vec{a}_t', \VEC(\mat{P}_t)' )'. \]

vector[] ssm_filter(vector[] y,
                    vector[] d, matrix[] Z, matrix[] H,
                    vector[] c, matrix[] T, matrix[] R, matrix[] Q,
                    vector a1, matrix P1) {

  // returned data
  vector[ssm_filter_size(dims(Z)[3], dims(Z)[2])] res[size(y)];
  int q;
  int n;
  int p;
  int m;

  // sizes
  n = size(y); // number of obs
  p = dims(Z)[2]; // obs size
  m = dims(Z)[3]; // number of states
  q = dims(Q)[2]; // number of state disturbances

  //print("Sizes: n = ", m, ", p = ", n, ", m = ", m, ", q = ", q);
  {
    // system matrices for current iteration
    vector[p] d_t;
    matrix[p, m] Z_t;
    matrix[p, p] H_t;
    vector[m] c_t;
    matrix[m, m] T_t;
    matrix[m, q] R_t;
    matrix[q, q] Q_t;
    matrix[m, m] RQR;
    // result matricees for each iteration
    vector[m] a;
    matrix[m, m] P;
    vector[p] v;
    matrix[p, p] Finv;
    matrix[m, p] K;
    real ll;
    int idx[6, 3];

    idx = ssm_filter_idx(m, p);

    d_t = d[1];
    Z_t = Z[1];
    H_t = H[1];
    c_t = c[1];
    T_t = T[1];
    R_t = R[1];
    Q_t = Q[1];
    RQR = quad_form(Q_t, R_t);
    a = a1;
    P = P1;
    for (t in 1:n) {
      if (t > 1) {
        if (size(d) > 1) {
          d_t = d[t];
        }
        if (size(Z) > 1) {
          Z_t = Z[t];
        }
        if (size(H) > 1) {
          H_t = H[t];
        }
        if (size(c) > 1) {
          c_t = c[t];
        }
        if (size(T) > 1) {
          T_t = T[t];
        }
        if (size(R) > 1) {
          R_t = R[t];
        }
        if (size(Q) > 1) {
          Q_t = Q[t];
        }
        if (size(R) > 1 && size(Q) > 1) {
          RQR = quad_form(Q_t, R_t);
        }
      }
      // updating
      v = ssm_filter_update_v(y[t], a, d_t, Z_t);
      Finv = ssm_filter_update_Finv(P, Z_t, H_t);
      K = ssm_filter_update_K(P, T_t, Z_t, Finv);
      ll = ssm_filter_update_ll(v, Finv);
      // saving
      res[t, 1] = ll;
      res[t, idx[2, 2]:idx[2, 3]] = v;
      res[t, idx[3, 2]:idx[3, 3]] = symmat_to_vector(Finv);
      res[t, idx[4, 2]:idx[4, 3]] = to_vector(K);
      res[t, idx[5, 2]:idx[5, 3]] = a;
      res[t, idx[6, 2]:idx[6, 3]] = symmat_to_vector(P);
      // predict a_{t + 1}, P_{t + 1}
      if (t < n) {
        a = ssm_filter_update_a(a, c_t, T_t, v, K);
        P = ssm_filter_update_P(P, Z_t, T_t, RQR, K);
      }
    }
  }
  return res;
}

4.3.10 ssm_filter_states

Parameters:

• int m Number of states

Return Value: int The size of the vector

Length of the vectors returned by ssm_filter_states

int ssm_filter_states_size(int m) {
  int sz;
  sz = m + symmat_size(m);
  return sz;
}

4.3.11 ssm_filter_states_get_a

Parameters:

• vector x A vector returned by ssm_filter_states
• int m Number of states

Return Value: matrix An \(m \times 1\) vector with the filtered expected value of the state, \(\vec{a}_{t|t} = \E(\vec{\alpha}_t | \vec{y}_{1:t})\).

Extract \(a_{t|t}\) from the results of ssm_filter_states

vector ssm_filter_states_get_a(vector x, int m) {
  vector[m] a;
  a = x[ :m];
  return a;
}

4.3.12 ssm_filter_states_get_P

Parameters:

• vector x A vector returned by ssm_filter_states
• int m Number of states

Return Value: matrix An \(m \times m\) matrix with the filtered variance of the state, \(\mat{P}_{t|t} = \Var(\vec{\alpha}_t | \vec{y}_{1:t})\).

Extract \(P_{t|t}\) from the results of ssm_filter_states

matrix ssm_filter_states_get_P(vector x, int m) {
  matrix[m, m] P;
  P = vector_to_symmat(x[(m + 1): ], m);
  return P;
}

4.3.13 ssm_filter_states

Parameters:

• vector[] filter Results from ssm_filter
• matrix[] Z Design matrix, \(\mat{Z}_t\). An array of \(p \times m\) matrices.

Return Value: Array of size \(n\) of vectors.

Calculate filtered expected values and variances of the states

The filtering function ssm_filter returns the mean and variance of the predicted states, \(\vec{a}_t = \E(\vec{\alpha}_t | \vec{y}_{1:(t -1)})\) and \(\mat{P}_t = \Var(\vec{\alpha}_t | \vec{y}_{1:(t -1)})\).

The vectors returned by ssm_filter_states are of length \(m + m ^ 2\), with \[ \vec{v}_t = (\vec{a}_{t|t}', \VEC(\vec{P}_{t|t})' )' \] Use the functions ssm_filter_states_get_a and ssm_filter_states_get_P to extract elements from the results.

For Z the array can have a size of 1, if it is not time-varying, or a size of \(n - 1\) if it is time varying.

vector[] ssm_filter_states(vector[] filter, matrix[] Z) {
  vector[ssm_filter_states_size(dims(Z)[3])] res[size(filter)];
  int n;
  int m;
  int p;
  n = size(filter);
  m = dims(Z)[3];
  p = dims(Z)[2];
  {
    // system matrices for current iteration
    matrix[p, m] Z_t;
    // filter matrices
    vector[m] aa; // filtered values of the state, a_{t|t}
    matrix[m, m] PP; // filtered values of the variance of the state, P_{t|t}
    vector[p] v;
    matrix[p, p] Finv;
    vector[m] a;
    matrix[m, m] P;

    Z_t = Z[1];
    for (t in 1:n) {
      if (t > 1) {
        if (size(Z) > 1) {
          Z_t = Z[t];
        }
      }
      // extract values from the filter
      v = ssm_filter_get_v(filter[t], m, p);
      Finv = ssm_filter_get_Finv(filter[t], m, p);
      a = ssm_filter_get_a(filter[t], m, p);
      P = ssm_filter_get_P(filter[t], m, p);
      // calcualte filtered values
      aa = a + P * Z_t ' * Finv * v;
      PP = to_symmetric_matrix(P - P * quad_form(Finv, Z_t) * P);
      // saving
      res[t, :m] = aa;
      res[t, (m + 1): ] = symmat_to_vector(PP);
    }
  }
  return res;
}

4.4 Log-likelihood

4.4.1 ssm_lpdf

Parameters:

• vector[] y Observations, \(\vec{y}_t\). An array of size \(n\) of \(p \times 1\) vectors.
• vector[] d Observation intercept, \(\vec{d}_t\). An array of \(p \times 1\) vectors.
• matrix[] Z Design matrix, \(\mat{Z}_t\). An array of \(p \times m\) matrices.
• matrix[] H Observation covariance matrix, \(\mat{H}_t\). An array of \(p \times p\) matrices.
• vector[] c State intercept, \(\vec{c}_t\). An array of \(m \times 1\) vectors.
• matrix[] T Transition matrix, \(\mat{T}_t\). An array of \(m \times m\) matrices.
• matrix[] R State covariance selection matrix, \(\mat{R} _t\). An array of \(p \times q\) matrices.
• matrix[] Q State covariance matrix, \(\mat{Q}_t\). An array of \(q \times q\) matrices.
• vector a1 Expected value of the intial state, \(a_1 = \E(\alpha_1)\). An \(m \times 1\) matrix.
• matrix P1 Variance of the initial state, \(P_1 = \Var(\alpha_1)\). An \(m \times m\) matrix.

Return Value: real The log-likelihood, \(p(\vec{y}_{1:n} | \vec{d}, \mat{Z}, \mat{H}, \vec{c}, \mat{T}, \mat{R}, \mat{Q})\), marginalized over the latent states.

Log-likelihood of a Linear Gaussian State Space Model

For d, Z, H, c, T, R, Q the array can have a size of 1, if it is not time-varying, or a size of \(n\) (for d, Z, H) or \(n - 1\) (for c, T, R, Q) if it is time varying.

The log-likelihood of a linear Gaussian state space model is, If the the system matrices and initial conditions are known, the log likelihood is \[ \begin{aligned}[t] \log L(\mat{Y}_n) &= \log p(\vec{y}_1, \dots, \vec{y}_n) = \sum_{t = 1}^n \log p(\vec{y}_t | \mat{Y}_{t - 1}) \\ &= - \frac{np}{2} \log 2 \pi - \frac{1}{2} \sum_{t = 1}^n \left( \log \left| \mat{F}_t \right| + \vec{v}\T \mat{F}_t^{-1} \vec{v}_t \right) \end{aligned} , \] where \(\mat{F}_t\) and \(\mat{V}_t\) come from a forward pass of the Kalman filter.

real ssm_lpdf(vector[] y,
               vector[] d, matrix[] Z, matrix[] H,
               vector[] c, matrix[] T, matrix[] R, matrix[] Q,
               vector a1, matrix P1) {
  real ll;
  int n;
  int m;
  int p;
  int q;
  n = size(y); // number of obs
  m = dims(Z)[2];
  p = dims(Z)[3];
  q = dims(Q)[2];
  {
    // system matrices for current iteration
    vector[p] d_t;
    matrix[p, m] Z_t;
    matrix[p, p] H_t;
    vector[m] c_t;
    matrix[m, m] T_t;
    matrix[m, q] R_t;
    matrix[q, q] Q_t;
    matrix[m, m] RQR;
    // result matricees for each iteration
    vector[n] ll_obs;
    vector[m] a;
    matrix[m, m] P;
    vector[p] v;
    matrix[p, p] Finv;
    matrix[m, p] K;

    d_t = d[1];
    Z_t = Z[1];
    H_t = H[1];
    c_t = c[1];
    T_t = T[1];
    R_t = R[1];
    Q_t = Q[1];
    RQR = quad_form(Q_t, R_t);

    a = a1;
    P = P1;
    for (t in 1:n) {
      if (t > 1) {
        if (size(d) > 1) {
          d_t = d[t];
        }
        if (size(Z) > 1) {
          Z_t = Z[t];
        }
        if (size(H) > 1) {
          H_t = H[t];
        }
        if (size(c) > 1) {
          c_t = c[t];
        }
        if (size(T) > 1) {
          T_t = T[t];
        }
        if (size(R) > 1) {
          R_t = R[t];
        }
        if (size(Q) > 1) {
          Q_t = Q[t];
        }
        if (size(R) > 1 && size(Q) > 1) {
          RQR = quad_form(Q_t, R_t);
        }
      }
      v = ssm_filter_update_v(y[t], a, d_t, Z_t);
      Finv = ssm_filter_update_Finv(P, Z_t, H_t);
      K = ssm_filter_update_K(P, Z_t, T_t, Finv);
      ll_obs[t] = ssm_filter_update_ll(v, Finv);
      // don't save a, P for last iteration
      if (t < n) {
        a = ssm_filter_update_a(a, c_t, T_t, v, K);
        P = ssm_filter_update_P(P, Z_t, T_t, RQR, K);
      }
    }
    ll = sum(ll_obs);
  }
  return ll;
}

4.5 Time-Invariant Kalman Filter

4.5.1 ssm_check_matrix_equal

Parameters:

• matrix A An \(m \times n\) matrix.
• matrix B An \(m \times n\) matrix.
• real The relative tolerance for convergence.

Return Value: int If converged, then 1, else 0.

Check if two matrices are approximately equal

The matrices \(A\) and \(B\) are considered approximately equal if \[ \max(A - B) / \max(A) < \epsilon, \] where \(\epsilon\) is the tolerance.

int ssm_check_matrix_equal(matrix A, matrix B, real tol) {
  real eps;
  eps = max(to_vector(A - B)) / max(to_vector(A));
  if (eps < tol) {
    return 1;
  } else {
    return 0;
  }
}

4.5.2 ssm_constant_lpdf

Parameters:

• vector[] y Observations, \(\vec{y}_t\). An array of size \(n\) of \(p \times 1\) vectors.
• vector d Observation intercept, \(\vec{d}_t\). An array of \(p \times 1\) vectors.
• matrix Z Design matrix, \(\mat{Z}_t\). An array of \(p \times m\) matrices.
• matrix H Observation covariance matrix, \(\mat{H}_t\). An array of \(p \times p\) matrices.
• vector c State intercept, \(\vec{c}_t\). An array of \(m \times 1\) vectors.
• matrix T Transition matrix, \(\mat{T}_t\). An array of \(m \times m\) matrices.
• matrix R State covariance selection matrix, \(\mat{R} _t\). An array of \(p \times q\) matrices.
• matrix Q State covariance matrix, \(\mat{Q}_t\). An array of \(q \times q\) matrices.
• vector a1 Expected value of the intial state, \(a_1 = \E(\alpha_1)\). An \(m \times 1\) matrix.
• matrix P1 Variance of the initial state, \(P_1 = \Var(\alpha_1)\). An \(m \times m\) matrix.

Return Value: real The log-likelihood, \(p(\vec{y}_{1:n} | \vec{d}, \mat{Z}, \mat{H}, \vec{c}, \mat{T}, \mat{R}, \mat{Q})\), marginalized over the latent states.

Log-likelihood of a Time-Invariant Linear Gaussian State Space Model

Unlike ssm_filter, this function requires the system matrices (d, Z, H, c, T, R, Q) to all be time invariant (constant). When the state space model is time-invariant, then the Kalman recursion for \(\mat{P}_t\) converges. This function takes advantage of this feature and stops updating \(\mat{P}_t\) after it converges to a steady state.

real ssm_constant_lpdf(vector[] y,
                      vector d, matrix Z, matrix H,
                      vector c, matrix T, matrix R, matrix Q,
                      vector a1, matrix P1) {
  real ll;
  int n;
  int m;
  int p;

  n = size(y); // number of obs
  m = cols(Z);
  p = rows(Z);
  {
    vector[n] ll_obs;
    vector[m] a;
    matrix[m, m] P;
    vector[p] v;
    matrix[p, p] Finv;
    matrix[m, p] K;
    matrix[m, m] RQR;
    // indicator for if the filter has converged
    // This only works for time-invariant state space models
    int converged;
    matrix[m, m] P_old;
    real tol;
    converged = 0;
    tol = 1e-7;

    RQR = quad_form(Q, R);
    a = a1;
    P = P1;
    for (t in 1:n) {
      v = ssm_filter_update_v(y[t], a, d, Z);
      if (converged < 1) {
        Finv = ssm_filter_update_Finv(P, Z, H);
        K = ssm_filter_update_K(P, Z, T, Finv);
      }
      ll_obs[t] = ssm_filter_update_ll(v, Finv);
      // don't save a, P for last iteration
      if (t < n) {
        a = ssm_filter_update_a(a, c, T, v, K);
        // check for convergence
        // should only check for convergence if there are no missing values
        if (converged < 1) {
          P_old = P;
          P = ssm_filter_update_P(P, Z, T, RQR, K);
          converged = ssm_check_matrix_equal(P, P_old, tol);
        }
      }
    }
    ll = sum(ll_obs);
  }
  return ll;
}

4.6 Common Smoother Functions

4.6.1 ssm_smooth_update_r

Parameters:

• vector r An \(m \times 1\) vector with \(\vec{r}_{t-1}\)
• matrix Z A \(p \times m\) vector with \(\mat{Z}_{t}\)
• vector v A \(p \times 1\) vector of the forecast errors, \(\vec{v}_t\).
• matrix Finv A \(p \times p\) matrix of the forecast precision, \(\mat{F}^{-1}_t\).
• matrix L An \(m \times m\) matrix with \(\mat{L}_t\).

Return Value: matrix An \(m \times 1\) vector with \(\vec{r}_t\).

Update \(\vec{r}_t\) in smoothing recursions

In smoothing recursions, the vector \(\vec{r}_t\) is updated with, \[ \vec{r}_{t - 1} = \mat{Z}' \mat{F}^{-1}_t \vec{v}_t + \mat{L}' \vec{r}_{t} . \]

See (Durbin and Koopman 2012, 91)

vector ssm_smooth_update_r(vector r, matrix Z, vector v, matrix Finv,
                           matrix L) {
  vector[num_elements(r)] r_new;
  r_new = Z ' * Finv * v + L ' * r;
  return r_new;
}

4.6.2 ssm_smooth_update_N

Parameters:

• vector N An \(m \times 1\) vector with \(\vec{N}_{t-1}\)
• matrix Z A \(p \times m\) vector with \(\mat{Z}_{t}\)
• matrix Finv A \(p \times p\) matrix of the forecast precision, \(\mat{F}^{-1}_t\).
• matrix L An \(m \times m\) matrix with \(\mat{L}_t\).

Return Value: matrix An \(m \times m\) matrix with \(\vec{N}_t\).

Update \(\mat{N}_t\) in smoothing recursions

In smoothing recursions, the matrix \(\vec{N}_t\) is updated with, \[ \mat{N}_{t - 1} = \mat{Z}_t' \mat{F}^{-1}_t \mat{Z}_t + \mat{L}_t' \mat{N}_t \mat{L}_t . \]

See (Durbin and Koopman 2012, 91)

matrix ssm_smooth_update_N(matrix N, matrix Z, matrix Finv, matrix L) {
  matrix[rows(N), cols(N)] N_new;
  N_new = quad_form(Finv, Z) + quad_form(N, L);
  return N_new;
}

4.6.3 ssm_smooth_state_size

Parameters:

• int m The number of states.

Return Value: int The size of the vectors is \(m + m (m + 1) / 2\).

The number of elements in vectors returned by ssm_smooth_state

int ssm_smooth_state_size(int m) {
  int sz;
  sz = m + symmat_size(m);
  return sz;
}

4.6.4 ssm_smooth_state_get_mean

Parameters:

• vector x A vector returned by ssm_smooth_state
• int q The number of state disturbances, \(\vec{\eta}_t\).

Return Value: vector An \(m \times 1\) vector with \(\hat{\vec{\eta}}_t\).

Extract \(\hat{\vec{\alpha}}_t\) from vectors returned by ssm_smooth_state

vector ssm_smooth_state_get_mean(vector x, int m) {
  vector[m] alpha;
  alpha = x[ :m];
  return alpha;
}

4.6.5 ssm_smooth_state_get_var

Parameters:

• vector x A vector returned by ssm_smooth_state
• int m The number of states

Return Value: matrix An \(m \times m\) matrix with \(\mat{V}_t\).

Extract \(mat{V}_t\) from vectors returned by ssm_smooth_state

matrix ssm_smooth_state_get_var(vector x, int m) {
  matrix[m, m] V;
  V = vector_to_symmat(x[(m + 1): ], m);
  return V;
}

4.6.6 ssm_smooth_state

Parameters:

• vector[] filter Results of ssm_filter
• matrix[] Z Design matrix, \(\mat{Z}_t\). An array of \(p \times m\) matrices.
• matrix[] T Transition matrix, \(\mat{T}_t\). An array of \(m \times m\) matrices.

Return Value: vector[] An array of vectors constaining \(\hat{\vec{\alpha}}_t\) and \(\mat{V}_t = \Var(\vec{\alpha}_t | \vec{y}_{1:n})\).

The state smoother

This calculates the mean and variance of the states, \(\vec{\alpha}_t\), given the entire sequence, \(\vec{y}_{1:n}\).

in the format described below.

For Z and T the array can have a size of 1, if it is not time-varying, or a size of \(n\) (for Z) or \(n - 1\) (for T) if it is time varying.

The vectors returned by this function have \(m + m ^ 2\) elements in this format, \[ (\hat{\vec{\alpha}}_t', \VEC(\mat{V}_t)' )'. \] Use the ssm_smooth_state_get_mean and ssm_smooth_state_get_var to extract components from the returned vectors.

value	length	start	end
\(\hat{\vec{\alpha}}_t\)	\(m\)	\(1\)	\(m\)
\(\mat{V}_t\)	\(m (m + 1) / 2\)	\(m + 1\)	\(m + m (m + 1) / 2\)

See Durbin and Koopman (2012), Eq 4.44 and eq 4.69.

vector[] ssm_smooth_state(vector[] filter, matrix[] Z, matrix[] T) {
  vector[ssm_smooth_state_size(dims(Z)[3])] res[size(filter)];
  int n;
  int m;
  int p;
  n = size(filter);
  m = dims(Z)[3];
  p = dims(Z)[2];
  {
    // system matrices for current iteration
    matrix[p, m] Z_t;
    matrix[m, m] T_t;
    // smoother results
    vector[m] r;
    matrix[m, m] N;
    matrix[m, m] L;
    vector[m] alpha;
    matrix[m, m] V;
    // filter results
    vector[p] v;
    matrix[m, p] K;
    matrix[p, p] Finv;
    vector[m] a;
    matrix[m, m] P;

    if (size(Z) == 1) {
      Z_t = Z[1];
    }
    if (size(T) == 1) {
      T_t = T[1];
    }
    // initialize smoother
    // r and N go from n, n - 1, ..., 1, 0.
    // r_n and N_n
    r = rep_vector(0.0, m);
    N = rep_matrix(0.0, m, m);
    // move backwards in time: t, ..., 1
    for (i in 0:(n - 1)) {
      int t;
      t = n - i;
      // set time-varying system matrices
      if (size(Z) > 1) {
        Z_t = Z[t];
      }
      if (size(T) > 1) {
        T_t = T[t];
      }
      // get filtered values
      K = ssm_filter_get_K(filter[t], m, p);
      v = ssm_filter_get_v(filter[t], m, p);
      Finv = ssm_filter_get_Finv(filter[t], m, p);
      a = ssm_filter_get_a(filter[t], m, p);
      P = ssm_filter_get_P(filter[t], m, p);
      // updating
      // L_t
      L = ssm_filter_update_L(Z_t, T_t, K);
      // r_{t - 1} and N_{t - 1}
      r = ssm_smooth_update_r(r, Z_t, v, Finv, L);
      N = ssm_smooth_update_N(N, Z_t, Finv, L);
      // hat(alpha)_{t} and V_t which use r and N from (t - 1)
      alpha = a + P * r;
      V = to_symmetric_matrix(P - P * N * P);
      // saving
      res[t, :m] = alpha;
      res[t, (m + 1): ] = symmat_to_vector(V);
    }
  }
  return res;
}

4.6.7 ssm_smooth_eps_size

Parameters:

• int p The length of the observation vectors, \(\vec{y}_t\).

Return Value: int The size of the vectors is \(p + p (p + 1) / 2\).

The size of the vectors returned by ssm_smooth_eps

int ssm_smooth_eps_size(int p) {
  int sz;
  sz = p + symmat_size(p);
  return sz;
}

4.6.8 ssm_smooth_eps_get_mean

Parameters:

• x A vector from the results of ssm_smooth_eps.
• int p The length of the observation vectors, \(\vec{y}_t\).

Return Value: vector A \(p \times 1\) vector with \(\hat{\vec{\varepsilon}}_t\).

Extract \(\hat{\vec{\varepsilon}}_t\) from vectors returned by ssm_smooth_eps

vector ssm_smooth_eps_get_mean(vector x, int p) {
  vector[p] eps;
  eps = x[ :p];
  return eps;
}

4.6.9 ssm_smooth_eps_get_var

Parameters:

• vector x A vector returned by ssm_smooth_eps
• int p The length of the observation vectors, \(\vec{y}_t\).

Return Value: matrix A \(p \times p\) matrix with \(\Var(\vec{\varepsilon}_t | \vec{y}_{1:n})\)

Extract \(\Var(\varepsilon_t|\vec{y}_{1:n})\) from vectors returned by ssm_smooth_eps

matrix ssm_smooth_eps_get_var(vector x, int p) {
  matrix[p, p] eps_var;
  eps_var = vector_to_symmat(x[(p + 1): ], p);
  return eps_var;
}

4.6.10 ssm_smooth_eps

Parameters:

• vector[] filter Results of ssm_filter
• matrix[] Z Design matrix, \(\mat{Z}_t\). An array of \(p \times m\) matrices.
• matrix[] H Observation covariance matrix, \(\mat{H}_t\). An array of \(p \times p\) matrices.
• matrix[] T Transition matrix, \(\mat{T}_t\). An array of \(m \times m\) matrices.

Return Value: vector[] An array of vectors constaining \(\hat{\vec{\varepsilon}}_t\) and \(\Var(\vec{\varepsilon}_t | \vec{y}_{1:n})\) in the format described below.

The observation disturbance smoother

This calculates the mean and variance of the observation disturbances, \(\vec{\varepsilon}_t\), given the entire sequence, \(\vec{y}_{1:n}\).

For Z,H, T, the array can have a size of 1, if it is not time-varying, or a size of \(n\) (for Z, H) or \(n - 1\) (for T), if it is time varying.

The vectors returned by this function have \(p + p (p + 1) / 2\) elements in this format, \[ (\hat{\vec{\varepsilon}}_t', \VEC(\Var(\vec{\varepsilon}_t | \vec{y}_{1:n}))' )' \]

value	length	start	end
\(\hat{\vec{\varepsilon}}_t\)	\(p\)	\(1\)	\(p\)
\(\Var(\vec{\varepsilon}_t | \vec{y}_{1:n})\)	\(p (p + 1) / 2\)	\(p + 1\)	\(p + p (p + 1) / 2\)

See (Durbin and Koopman 2012, Sec 4.5.3 (eq 4.69))

vector[] ssm_smooth_eps(vector[] filter, matrix[] Z, matrix[] H, matrix[] T) {
  vector[ssm_smooth_eps_size(dims(Z)[2])] res[size(filter)];
  int n;
  int m;
  int p;
  n = size(filter);
  m = dims(Z)[3];
  p = dims(Z)[2];
  {
    // smoother values
    vector[m] r;
    matrix[m, m] N;
    matrix[m, m] L;
    vector[p] eps;
    matrix[p, p] var_eps;
    // filter results
    vector[p] v;
    matrix[m, p] K;
    matrix[p, p] Finv;
    // system matrices
    matrix[p, m] Z_t;
    matrix[p, p] H_t;
    matrix[m, m] T_t;

    // set matrices if time-invariant
    if (size(Z) == 1) {
      Z_t = Z[1];
    }
    if (size(H) == 1) {
      H_t = H[1];
    }
    if (size(T) == 1) {
      T_t = T[1];
    }
    // initialize smoother
    // r and N go from n, n - 1, ..., 1, 0.
    // r_n and N_n
    r = rep_vector(0.0, m);
    N = rep_matrix(0.0, m, m);
    for (i in 1:n) {
      int t;
      // move backwards in time
      t = n - i + 1;
      // update time-varying system matrices
      if (size(Z) > 1) {
        Z_t = Z[t];
      }
      if (size(H) > 1) {
        H_t = H[t];
      }
      if (size(T) > 1) {
        T_t = T[t];
      }
      // get values from filter
      K = ssm_filter_get_K(filter[t], m, p);
      v = ssm_filter_get_v(filter[t], m, p);
      Finv = ssm_filter_get_Finv(filter[t], m, p);
      // updating
      L = ssm_filter_update_L(Z_t, T_t, K);
      // r_{t - 1} and N_{t - 1}
      r = ssm_smooth_update_r(r, Z_t, v, Finv, L);
      N = ssm_smooth_update_N(N, Z_t, Finv, L);
      // eps_t and V(eps_t|y)
      eps = H_t * (Finv * v - K ' * r);
      var_eps = to_symmetric_matrix(H_t - H_t * (Finv + quad_form(N, K)) * H_t);
      // saving
      res[t, :p] = eps;
      res[t, (p + 1): ] = symmat_to_vector(var_eps);
    }
  }
  return res;
}

4.6.11 ssm_smooth_eta

Parameters:

• int p The length of the observation vectors, \(\vec{y}_t\).

Return Value: int The size of the vectors is \(q + q (q + 1) / 2\).

The size of the vectors returned by ssm_smooth_eta

int ssm_smooth_eta_size(int q) {
  int sz;
  sz = q + symmat_size(q);
  return sz;
}

4.6.12 ssm_smooth_eta_get_mean

Parameters:

• vector x A vector returned by ssm_smooth_eta
• int q The number of state disturbances, \(\vec{\eta}_t\).

Return Value: vector A \(q \times 1\) vector with \(\hat{\vec{\eta}}_t\).

Extract \(\hat{\vec{\varepsilon}}_t\) from vectors returned by ssm_smooth_eta

vector ssm_smooth_eta_get_mean(vector x, int q) {
  vector[q] eta;
  eta = x[ :q];
  return eta;
}

4.6.13 ssm_smooth_eta_get_var

Parameters:

• vector x A vector returned by ssm_smooth_eta
• int q The number of state disturbances, \(\vec{\eta}_t\).

Return Value: matrix A \(q \times q\) matrix with \(\Var(\vec{\eta}_t | \vec{y}_{1:n})\).

Extract \(\Var(\eta_t|\vec{y}_{1:n})\) from vectors returned by ssm_smooth_eta

matrix ssm_smooth_eta_get_var(vector x, int q) {
  matrix[q, q] eta_var;
  eta_var = vector_to_symmat(x[(q + 1): ], q);
  return eta_var;
}

4.6.14 ssm_smooth_eta

Parameters:

• vector[] filter Results of ssm_filter
• matrix[] Z Design matrix, \(\mat{Z}_t\). An array of \(p \times m\) matrices.
• matrix[] T Transition matrix, \(\mat{T}_t\). An array of \(m \times m\) matrices.
• matrix[] R State covariance selection matrix, \(\mat{R} _t\). An array of \(p \times q\) matrices.
• matrix[] Q State covariance matrix, \(\mat{Q}_t\). An array of \(q \times q\) matrices.

Return Value: vector[] An array of vectors constaining \(\hat{\vec{\eta}}_t\) and \(\Var(\vec{\eta}_t | \vec{y}_{1:n})\) in the format described below.

The state disturbance smoother

This calculates the mean and variance of the observation disturbances, \(\vec{\eta}_t\), given the entire sequence, \(\vec{y}_{1:n}\).

For Z, T, R, Q the array can have a size of 1, if it is not time-varying, or a size of \(n\) (for Z) or \(n - 1\) (for T, R, Q) if it is time varying.

The vectors returned by this function have \(q + q (q + 1) / 2\) elements in this format, \[ (\hat{\vec{\eta}}_t', \VEC(\Var(\vec{\eta}_t | \vec{y}_{1:n}))' ). \] Use the ssm_smooth_eta_get_mean and ssm_smooth_eta_get_var to extract components from the returned vectors.

value	length	start	end
\(\hat{\vec{\eta}}_t\)	\(q\)	\(1\)	\(q\)
\(\Var(\vec{\eta}_t | \vec{y}_{1:n})\)	\(q (q + 1) / 2\)	\(q + 1\)	\(q + q (q + 1) / 2\)

See (Durbin and Koopman 2012, Sec 4.5.3 (eq 4.69))

vector[] ssm_smooth_eta(vector[] filter,
                        matrix[] Z, matrix[] T,
                        matrix[] R, matrix[] Q) {
  vector[ssm_smooth_eta_size(dims(Q)[2])] res[size(filter)];
  int n;
  int m;
  int p;
  int q;
  n = size(filter);
  m = dims(Z)[3];
  p = dims(Z)[2];
  q = dims(Q)[2];
  {
    // smoother matrices
    vector[m] r;
    matrix[m, m] N;
    matrix[m, m] L;
    vector[q] eta;
    matrix[q, q] var_eta;
    // system matrices
    matrix[p, m] Z_t;
    matrix[m, m] T_t;
    matrix[m, q] R_t;
    matrix[q, q] Q_t;
    // filter matrices
    vector[p] v;
    matrix[m, p] K;
    matrix[p, p] Finv;

    // set time-invariant matrices
    if (size(Z) == 1) {
      Z_t = Z[1];
    }
    if (size(T) == 1) {
      T_t = T[1];
    }
    if (size(R) == 1) {
      R_t = R[1];
    }
    if (size(Q) == 1) {
      Q_t = Q[1];
    }
    // initialize smoother
    r = rep_vector(0.0, m);
    N = rep_matrix(0.0, m, m);
    for (i in 0:(n - 1)) {
      int t;
      // move backwards in time
      t = n - i;
      // update time-varying system matrices
      if (size(Z) > 1) {
        Z_t = Z[t];
      }
      if (size(T) > 1) {
        T_t = T[t];
      }
      if (size(R) > 1) {
        R_t = R[t];
      }
      if (size(Q) > 1) {
        Q_t = Q[t];
      }
      // get values from filter
      K = ssm_filter_get_K(filter[t], m, p);
      v = ssm_filter_get_v(filter[t], m, p);
      Finv = ssm_filter_get_Finv(filter[t], m, p);
      // update smoother
      L = ssm_filter_update_L(Z_t, T_t, K);
      r = ssm_smooth_update_r(r, Z_t, v, Finv, L);
      N = ssm_smooth_update_N(N, Z_t, Finv, L);
      eta = Q_t * R_t ' * r;
      var_eta = to_symmetric_matrix(Q_t - Q_t * quad_form(N, R_t) * Q_t);
      // saving
      res[t, :q] = eta;
      res[t, (q + 1): ] = symmat_to_vector(var_eta);
    }
  }
  return res;
}

4.6.15 ssm_smooth_faststate

Parameters:

• vector[] filter The results of ssm_filter
• matrix[] Z Design matrix, \(\mat{Z}_t\). An array of \(p \times m\) matrices.
• vector[] c State intercept, \(\vec{c}_t\). An array of \(m \times 1\) vectors.
• matrix[] T Transition matrix, \(\mat{T}_t\). An array of \(m \times m\) matrices.
• matrix[] R State covariance selection matrix, \(\mat{R} _t\). An array of \(p \times q\) matrices.
• matrix[] Q State covariance matrix, \(\mat{Q}_t\). An array of \(q \times q\) matrices.

Return Value: vector[] An array of size \(n\) of \(m \times 1\) vectors containing \(\hat{\vec{\alpha}}_t\).

The fast state smoother

The fast state smoother calculates \(\hat{\vec{\alpha}}_t = \E(\vec{\alpha}_t | \vec{y}_{1:n})\). \[ \hat{\vec{\alpha}}_{t + 1} = \mat{T}_t \hat{\vec{\alpha}}_{t} + \mat{R}_t \mat{Q}_t \mat{R}'_t \vec{r}_t , \] where \(r_t\) is calcualted from the state disturbance smoother. The smoother is initialized at \(t = 1\) with \(\hat{\vec{\alpha}}_t = \vec{a}_1 + \mat{P}_1 \vec{r}_0\).

Unlike the normal state smoother, it does not calculate the variances of the smoothed state.

For Z, c, T, R, Q the array can have a size of 1, if it is not time-varying, or a size of \(n\) (for Z) or \(n - 1\) (for c, T, R, Q) if it is time varying.

See (Durbin and Koopman 2012, Sec 4.5.3 (eq 4.69))

vector[] ssm_smooth_faststate(vector[] filter,
                              vector[] c, matrix[] Z, matrix[] T,
                              matrix[] R, matrix[] Q) {
  vector[dims(Z)[3]] alpha[size(filter)];
  int n;
  int m;
  int p;
  int q;
  n = size(filter);
  m = dims(Z)[3];
  p = dims(Z)[2];
  q = dims(Q)[2];
  {
    // smoother matrices
    vector[m] r[n + 1];
    matrix[m, m] L;
    vector[m] a1;
    matrix[m, m] P1;
    // filter matrices
    vector[p] v;
    matrix[m, p] K;
    matrix[p, p] Finv;
    // system matrices
    matrix[p, m] Z_t;
    vector[m] c_t;
    matrix[m, m] T_t;
    matrix[p, q] R_t;
    matrix[q, q] Q_t;
    matrix[m, m] RQR;
    // set time-invariant matrices
    if (size(c) == 1) {
      c_t = c[1];
    }
    if (size(Z) == 1) {
      Z_t = Z[1];
    }
    if (size(T) == 1) {
      T_t = T[1];
    }
    if (size(R) == 1) {
      R_t = R[1];
    }
    if (size(Q) == 1) {
      Q_t = Q[1];
    }
    if (size(Q) == 1 && size(R) == 1) {
      RQR = quad_form(Q[1], R[1]');
    }
    // find smoothed state disturbances
    // Since I don't need to calculate the
    // variances of the smoothed disturbances,
    // I reimplement the state distrurbance smoother here
    // removing extraneous parts.
    // r goes from t = n, ..., 1, 0.
    // r_n
    r[n + 1] = rep_vector(0.0, m);
    for (i in 0:(n - 1)) {
      int t;
      // move backwards in time
      t = n - i;
      // update time varying system matrices
      if (size(Z) > 1) {
        Z_t = Z[t];
      }
      if (size(T) > 1) {
        T_t = T[t];
      }
      // get filter values
      K = ssm_filter_get_K(filter[t], m, p);
      v = ssm_filter_get_v(filter[t], m, p);
      Finv = ssm_filter_get_Finv(filter[t], m, p);
      // updating smoother
      L = ssm_filter_update_L(Z_t, T_t, K);
      // r_{t - 1}
      r[t] = ssm_smooth_update_r(r[t + 1], Z_t, v, Finv, L);
    }
    // calculate smoothed states
    a1 = ssm_filter_get_a(filter[1], m, p);
    P1 = ssm_filter_get_P(filter[1], m, p);
    // r[1] = r_0
    alpha[1] = a1 + P1 * r[1];
    // 1:(n - 1) -> \alpha_{2}:\alpha_{n}
    for (t in 1:(n - 1)) {
      if (size(c) > 1) {
        c_t = c[t];
      }
      if (size(T) > 1) {
        T_t = T[t];
      }
      if (size(Q) > 1) {
        Q_t = Q[t];
      }
      if (size(R) > 1) {
        R_t = R[t];
      }
      if (size(Q) > 1 || size(R) > 1) {
        RQR = quad_form(Q_t, R_t');
      }
      // `r[t + 1]` = $r_{t}$
      // alpha_{t + 1} = c_t + T_t * \alpha_t + R_t Q_t R'_t r_t
      alpha[t + 1] = c_t + T_t * alpha[t] + RQR * r[t + 1];
    }
  }
  return alpha;
}

4.7 Simulators and Smoothing Simulators

4.7.1 ssm_sim_idx

Parameters:

• int m The number of states
• int p The length of the observation vector
• int q The number of state disturbances

Return Value: int[,] A 4 x 3 array of integers

Indexes of each component of ssm_sim_rng results.

The returned array has columns (length, start location, and end location) for rows: \(\vec{y}_t\), \(\vec{\alpha}_t\), \(\vec{\varepsilon}_t\), and \(\vec{\eta}_t\) in the results of ssm_sim_rng.

int[,] ssm_sim_idx(int m, int p, int q) {
  int sz[4, 3];
  // y
  sz[1, 1] = p;
  // a
  sz[2, 1] = m;
  // eps
  sz[3, 1] = p;
  // eta
  sz[4, 1] = q;
  // Fill in start and stop points
  sz[1, 2] = 1;
  sz[1, 3] = sz[1, 2] + sz[1, 1] - 1;
  for (i in 2:4) {
    sz[i, 2] = sz[i - 1, 3] + 1;
    sz[i, 3] = sz[i, 2] + sz[i, 1] - 1;
  }
  return sz;
}

4.7.2 ssm_sim_size

Parameters:

• int m The number of states
• int p The length of the observation vector
• int q The number of state disturbances

Return Value: int The number of elements

The number of elements in vectors returned by ssm_sim_rng results.

int ssm_sim_size(int m, int p, int q) {
  int sz;
  sz = ssm_sim_idx(m, p, q)[4, 3];
  return sz;
}

4.7.3 ssm_sim_get_y

Parameters:

• int m The number of states
• int p The length of the observation vector
• int q The number of state disturbances

Return Value: vector vector A \(p \times 1\) vector with \(\vec{y}_t\).

Extract \(\vec{y}_t\) from vectors returned by ssm_sim_rng.

vector ssm_sim_get_y(vector x, int m, int p, int q) {
  vector[m] y;
  int idx[4, 3];
  idx = ssm_sim_idx(m, p, q);
  y = x[idx[1, 2]:idx[1, 3]];
  return y;
}

4.7.4 ssm_sim_get_a

Parameters:

• int m The number of states
• int p The length of the observation vector
• int q The number of state disturbances

Return Value: vector A \(m \times 1\) vector with \(\vec{\alpha}_t\).

Extract \(\vec{\alpha}_t\) from vectors returne by ssm_sim_rng.

vector ssm_sim_get_a(vector x, int m, int p, int q) {
  vector[m] a;
  int idx[4, 3];
  idx = ssm_sim_idx(m, p, q);
  a = x[idx[2, 2]:idx[2, 3]];
  return a;
}

4.7.5 ssm_sim_get_eps

Parameters:

• int m The number of states
• int p The length of the observation vector
• int q The number of state disturbances

Return Value: vector vector A \(p \times 1\) vector with \(\vec{\varepsilon}_t\).

Extract \(\vec{\varepsilon}_t\) from vectors returne by ssm_sim_rng.

vector ssm_sim_get_eps(vector x, int m, int p, int q) {
  vector[m] eps;
  int idx[4, 3];
  idx = ssm_sim_idx(m, p, q);
  eps = x[idx[3, 2]:idx[3, 3]];
  return eps;
}

4.7.6 ssm_sim_get_eta

Parameters:

• int m The number of states
• int p The length of the observation vector
• int q The number of state disturbances

Return Value: vector vector A \(q \times 1\) vector with \(\vec{\eta}_t\).

Extract \(\vec{\eta}_t\) from vectors returne by ssm_sim_rng.

vector ssm_sim_get_eta(vector x, int m, int p, int q) {
  vector[m] eta;
  int idx[4, 3];
  idx = ssm_sim_idx(m, p, q);
  eta = x[idx[4, 2]:idx[4, 3]];
  return eta;
}

4.7.7 ssm_sim_rng

Parameters:

• vector[] y Observations, \(\vec{y}_t\). An array of size \(n\) of \(p \times 1\) vectors.
• vector[] d Observation intercept, \(\vec{d}_t\). An array of \(p \times 1\) vectors.
• matrix[] Z Design matrix, \(\mat{Z}_t\). An array of \(p \times m\) matrices.
• matrix[] H Observation covariance matrix, \(\mat{H}_t\). An array of \(p \times p\) matrices.
• vector[] c State intercept, \(\vec{c}_t\). An array of \(m \times 1\) vectors.
• matrix[] T Transition matrix, \(\mat{T}_t\). An array of \(m \times m\) matrices.
• matrix[] R State covariance selection matrix, \(\mat{R} _t\). An array of \(p \times q\) matrices.
• matrix[] Q State covariance matrix, \(\mat{Q}_t\). An array of \(q \times q\) matrices.
• vector a1 Expected value of the intial state, \(a_1 = \E(\alpha_1)\). An \(m \times 1\) matrix.
• matrix P1 Variance of the initial state, \(P_1 = \Var(\alpha_1)\). An \(m \times m\) matrix.

Return Value: Array of size \(n\) of vectors with Draw \(\vec{y}_t\), \(\vec{\alpha}_t\), \(\vec{\eta}_t\) and \(\vec{\varepsilon}_t\). See the description.

Simulate from a Linear Gaussian State Space model.

For d, Z, H, c, T, R, Q the array can have a size of 1, if it is not time-varying, or a size of \(n\) (for d, Z, H) or \(n - 1\) (for c, T, R, Q) if it is time varying.

Draw \(\vec{y}_t\), \(\vec{\alpha}_t\), \(\vec{\eta}_t\) and \(\vec{\varepsilon}_t\) from the state space model, \[ \begin{aligned}[t] \vec{y}_t &= \vec{d}_t + \mat{Z}_t \vec{\alpha}_t + \vec{\varepsilon}_t, & \vec{\varepsilon}_t & \sim N(0, \mat{H}_t), \\ \vec{\alpha}_{t + 1} &= \vec{c}_t + \mat{T}_t \vec{\alpha}_t + \mat{R}_t \vec{\eta}_t, & \vec{\eta}_t & \sim N(0, \mat{Q}_t), \\ && \vec{\alpha}_1 &\sim N(\vec{a}_1, \mat{P}_1) . \end{aligned} \]

The returned vectors are of length \(2 p + m + q\), in the format, \[ (\vec{y}_t', \vec{\alpha}_t', \vec{\varepsilon}_t', \vec{\eta}_t') . \] Note that \(\eta_n = \vec{0}_q\). Use the functions ssm_sim_get_y, ssm_sim_get_a, ssm_sim_get_eps, and ssm_sim_get_eta to extract values from the vector.

element	length	start	end
\(y_t\)	\(p\)	\(1\)	\(p\)
\(\alpha\)_t	\(m\)	\(p + 1\)	\(p + m\)
\(\varepsilon_t\)	\(p\)	\(p + m + 1\)	\(2 p + m\)
\(\eta_t\)	\(q\)	\(2 p + m + 1\)	\(2 p + m + q\)

It is preferrable to use ssm_sim_get_y, ssm_sim_get_a, ssm_sim_get_eps, and ssm_sim_get_eta to extract values from these vectors.

vector[] ssm_sim_rng(int n,
                    vector[] d, matrix[] Z, matrix[] H,
                    vector[] c, matrix[] T, matrix[] R, matrix[] Q,
                    vector a1, matrix P1) {
  vector[ssm_sim_size(dims(Z)[3], dims(Z)[2], dims(Q)[2])] ret[n];
  int p;
  int m;
  int q;
  p = dims(Z)[2];
  m = dims(Z)[3];
  q = dims(Q)[2];
  {
    // system matrices for current iteration
    vector[p] d_t;
    matrix[p, m] Z_t;
    matrix[p, p] H_t;
    vector[m] c_t;
    matrix[m, m] T_t;
    matrix[m, q] R_t;
    matrix[q, q] Q_t;
    matrix[m, m] RQR;
    // outputs
    vector[p] y;
    vector[p] eps;
    vector[m] a;
    vector[q] eta;
    // constants
    vector[p] zero_p;
    vector[q] zero_q;
    vector[m] zero_m;
    int idx[4, 3];

    d_t = d[1];
    Z_t = Z[1];
    H_t = H[1];
    c_t = c[1];
    T_t = T[1];
    R_t = R[1];
    Q_t = Q[1];

    idx = ssm_sim_idx(m, p, q);
    zero_p = rep_vector(0.0, p);
    zero_q = rep_vector(0.0, q);
    zero_m = rep_vector(0.0, m);
    a = multi_normal_rng(a1, P1);
    for (t in 1:n) {
      // set system matrices
      if (t > 1) {
        if (size(d) > 1) {
          d_t = d[t];
        }
        if (size(Z) > 1) {
          Z_t = Z[t];
        }
        if (size(H) > 1) {
          H_t = H[t];
        }
        // system matrices are n - 1 length
        if (t < n) {
          if (size(c) > 1) {
            c_t = c[t];
          }
          if (size(T) > 1) {
            T_t = T[t];
          }
          if (size(R) > 1) {
            R_t = R[t];
          }
          if (size(Q) > 1) {
            Q_t = Q[t];
          }
        }
      }
      // draw forecast error
      eps = multi_normal_rng(zero_p, H_t);
      // draw observed value
      y = d_t + Z_t * a + eps;
      // since eta_t is for alpha_{t + 1}, we don't
      // draw it for t == n
      if (t == n) {
        eta = zero_q;
      } else {
        eta = multi_normal_rng(zero_q, Q_t);
      }
      // save
      ret[t, idx[1, 2]:idx[1, 3]] = y;
      ret[t, idx[2, 2]:idx[2, 3]] = a;
      ret[t, idx[3, 2]:idx[3, 3]] = eps;
      ret[t, idx[4, 2]:idx[4, 3]] = eta;
      // a_{t + 1}
      if (t < n) {
        a = c_t + T_t * a + R_t * eta;
      }
    }
  }
  return ret;
}

4.8 Simulation Smoothers

4.8.1 ssm_simsmo_state_rng

Parameters:

• vector[] alpha An of size \(n\) of \(m \times 1\) vectors containing the smoothed expected values of the states, \(\E(\vec{\alpha}_{1:n} | \vec{y}_{1:n})\). These are returned by sim_smooth_faststates. If sim_smooth_state was used, then the expected values need to first be extracted using sim_smooth_state_get_mean.
• vector[] d Observation intercept, \(\vec{d}_t\). An array of \(p \times 1\) vectors.
• matrix[] Z Design matrix, \(\mat{Z}_t\). An array of \(p \times m\) matrices.
• matrix[] H Observation covariance matrix, \(\mat{H}_t\). An array of \(p \times p\) matrices.
• vector[] c State intercept, \(\vec{c}_t\). An array of \(m \times 1\) vectors.
• matrix[] T Transition matrix, \(\mat{T}_t\). An array of \(m \times m\) matrices.
• matrix[] R State covariance selection matrix, \(\mat{R} _t\). An array of \(p \times q\) matrices.
• matrix[] Q State covariance matrix, \(\mat{Q}_t\). An array of \(q \times q\) matrices.
• vector a1 Expected value of the intial state, \(a_1 = \E(\alpha_1)\). An \(m \times 1\) matrix.
• matrix P1 Variance of the initial state, \(P_1 = \Var(\alpha_1)\). An \(m \times m\) matrix.

Return Value: vector[] Array of size \(n\) of \(m \times 1\) vectors containing a single draw from \((\vec{\alpha}_{1:n} | \vec{y}_{1:n})\).

State simulation smoother

Draw samples from the posterior distribution of the states, \(\tilde{\vec{\alpha}}_{1:n} \sim p(\vec{\alpha}_{1:n} | \vec{y}_{1:n})\).

For d, Z, H, c, T, R, Q the array can have a size of 1, if it is not time-varying, or a size of \(n\) (for d, Z, H) or \(n - 1\) (for c, T, R, Q) if it is time varying.

This draws samples using mean-correction simulation smoother of (Durbin and Koopman 2002). See (Durbin and Koopman 2012, Sec 4.9).

vector[] ssm_simsmo_states_rng(vector[] alpha,
                      vector[] d, matrix[] Z, matrix[] H,
                      vector[] c, matrix[] T, matrix[] R, matrix[] Q,
                      vector a1, matrix P1) {
    vector[dims(Z)[2]] draws[size(alpha)];
    int n;
    int p;
    int m;
    int q;
    n = size(alpha);
    p = dims(Z)[2];
    m = dims(Z)[3];
    q = dims(Q)[2];
    {
      vector[ssm_filter_size(m, p)] filter[n];
      vector[ssm_sim_size(m, p, q)] sims[n];
      vector[p] y[n];
      vector[m] alpha_hat_plus[n];
      // simulate unconditional disturbances and observations
      sims = ssm_sim_rng(n, d, Z, H, c, T, R, Q, a1, P1);
      for (i in 1:n) {
        y[i] = ssm_sim_get_y(sims[i], m, p, q);
      }
      // filter with simulated y's
      filter = ssm_filter(y, d, Z, H, c, T, R, Q, a1, P1);
      // mean correct epsilon samples
      alpha_hat_plus = ssm_smooth_faststate(filter, c, Z, T, R, Q);
      for (i in 1:n) {
        draws[i] = (ssm_sim_get_a(sims[i], m, p, q)
                    - alpha_hat_plus[i]
                    + alpha[i]);
      }
    }
    return draws;
}

4.8.2 ssm_simsmo_eta_rng

Parameters:

• vector[] eta Values returned by sim_smooth_eta
• vector[] d Observation intercept, \(\vec{d}_t\). An array of \(p \times 1\) vectors.
• matrix[] Z Design matrix, \(\mat{Z}_t\). An array of \(p \times m\) matrices.
• matrix[] H Observation covariance matrix, \(\mat{H}_t\). An array of \(p \times p\) matrices.
• vector[] c State intercept, \(\vec{c}_t\). An array of \(m \times 1\) vectors.
• matrix[] T Transition matrix, \(\mat{T}_t\). An array of \(m \times m\) matrices.
• matrix[] R State covariance selection matrix, \(\mat{R} _t\). An array of \(p \times q\) matrices.
• matrix[] Q State covariance matrix, \(\mat{Q}_t\). An array of \(q \times q\) matrices.
• vector a1 Expected value of the intial state, \(a_1 = \E(\alpha_1)\). An \(m \times 1\) matrix.
• matrix P1 Variance of the initial state, \(P_1 = \Var(\alpha_1)\). An \(m \times m\) matrix.

Return Value: vector[] Array of size \(n\) of \(q \times 1\) vectors containing a single draw from \((\vec{\eta}_{1:n} | \vec{y}_{1:n})\).

State disturbance simulation smoother

Draw samples from the posterior distribution of the observation disturbances, \(\tilde{\vec{\eta}}_{1:n} \sim p(\vec{\eta}_{1:n} | \vec{y}_{1:n})\).

For d, Z, H, c, T, R, Q the array can have a size of 1, if it is not time-varying, or a size of \(n\) (for d, Z, H) or \(n - 1\) (for c, T, R, Q) if it is time varying.

This draws samples using mean-correction simulation smoother of (Durbin and Koopman 2002). See (Durbin and Koopman 2012, Sec 4.9).

vector[] ssm_simsmo_eta_rng(vector[] eta,
                            vector[] d, matrix[] Z, matrix[] H,
                            vector[] c, matrix[] T, matrix[] R, matrix[] Q,
                            vector a1, matrix P1) {
    vector[dims(Q)[2]] draws[size(eta)];
    int n;
    int p;
    int m;
    int q;
    n = size(eta);
    p = dims(Z)[2];
    m = dims(Z)[3];
    q = dims(Q)[2];
    {
      vector[ssm_filter_size(m, p)] filter[n];
      vector[p] y[n];
      vector[ssm_sim_size(m, p, q)] sims[n];
      vector[ssm_smooth_eta_size(q)] etahat_plus[n];
      // simulate unconditional disturbances and observations
      sims = ssm_sim_rng(n, d, Z, H, c, T, R, Q, a1, P1);
      for (i in 1:n) {
        y[i] = ssm_sim_get_y(sims[i], m, p, q);
      }
      // filter simulated y's
      filter = ssm_filter(y, d, Z, H, c, T, R, Q, a1, P1);
      // mean correct eta samples
      etahat_plus = ssm_smooth_eta(filter, Z, T, R, Q);
      for (i in 1:n) {
        draws[i] = (ssm_sim_get_eta(sims[i], m, p, q)
                                    - ssm_smooth_eta_get_mean(etahat_plus[i], q)
                                    + ssm_smooth_eta_get_mean(eta[i], q));
      }
    }
    return draws;
}

4.8.3 ssm_simsmo_eps_rng

Parameters:

• vector[] eps Values returned by sim_smooth_eps
• vector[] d Observation intercept, \(\vec{d}_t\). An array of \(p \times 1\) vectors.
• matrix[] Z Design matrix, \(\mat{Z}_t\). An array of \(p \times m\) matrices.
• matrix[] H Observation covariance matrix, \(\mat{H}_t\). An array of \(p \times p\) matrices.
• vector[] c State intercept, \(\vec{c}_t\). An array of \(m \times 1\) vectors.
• matrix[] T Transition matrix, \(\mat{T}_t\). An array of \(m \times m\) matrices.
• matrix[] R State covariance selection matrix, \(\mat{R} _t\). An array of \(p \times q\) matrices.
• matrix[] Q State covariance matrix, \(\mat{Q}_t\). An array of \(q \times q\) matrices.
• vector a1 Expected value of the intial state, \(a_1 = \E(\alpha_1)\). An \(m \times 1\) matrix.
• matrix P1 Variance of the initial state, \(P_1 = \Var(\alpha_1)\). An \(m \times m\) matrix.

Return Value: vector[] Array of size \(n\) of \(p \times 1\) vectors containing a single draw from \((\vec{\varepsilon}_{1:n} | \vec{y}_{1:n})\).

Observation disturbance simulation smoother

Draw samples from the posterior distribution of the observation disturbances, \(\tilde{\vec{\varepsilon}}_{1:n} \sim p(\vec{\varepsilon}_{1:n} | \vec{y}_{1:n})\).

For d, Z, H, c, T, R, Q the array can have a size of 1, if it is not time-varying, or a size of \(n\) (for d, Z, H) or \(n - 1\) (for c, T, R, Q) if it is time varying.

This draws samples using mean-correction simulation smoother of (Durbin and Koopman 2002). See (Durbin and Koopman 2012, Sec 4.9).

vector[] ssm_simsmo_eps_rng(vector[] eps,
                      vector[] d, matrix[] Z, matrix[] H,
                      vector[] c, matrix[] T, matrix[] R, matrix[] Q,
                      vector a1, matrix P1) {
    vector[dims(Z)[2]] draws[size(eps)];
    int n;
    int p;
    int m;
    int q;
    n = size(eps);
    p = dims(Z)[2];
    m = dims(Z)[3];
    q = dims(Q)[2];
    {
      vector[ssm_filter_size(m, p)] filter[n];
      vector[p] y[n];
      vector[ssm_sim_size(m, p, q)] sims[n];
      vector[ssm_smooth_eta_size(p)] epshat_plus[n];
      // simulate unconditional disturbances and observations
      sims = ssm_sim_rng(n, d, Z, H, c, T, R, Q, a1, P1);
      for (i in 1:n) {
        y[i] = ssm_sim_get_y(sims[i], m, p, q);
      }
      // filter simulated y's
      filter = ssm_filter(y, d, Z, H, c, T, R, Q, a1, P1);
      // mean correct epsilon samples
      epshat_plus = ssm_smooth_eps(filter, Z, H, T);
      for (i in 1:n) {
        draws[i] = (ssm_sim_get_eps(sims[i], m, p, q)
                    - ssm_smooth_eps_get_mean(epshat_plus[i], p)
                    + ssm_smooth_eps_get_mean(eps[i], p));
      }
    }
    return draws;
}

4.9 Stationary

4.9.1 pacf_to_acf

Parameters:

• vector x A vector of coefficients of a partial autocorrelation function

Return Value: vector A vector of coefficients of an Autocorrelation function

Partial Autocorrelations to Autocorrelations

vector pacf_to_acf(vector x) {
  matrix[num_elements(x), num_elements(x)] y;
  int n;
  n = num_elements(x);
  y = rep_matrix(0.0, n, n);
  for (k in 1:n) {
    for (i in 1:(k - 1)) {
      y[k, i] = y[k - 1, i] + x[k] * y[k - 1, k - i];
    }
    y[k, k] = x[k];
    print(y);
  }
  return -y[n] ';
}

4.9.2 constrain_stationary

Parameters:

• vector x An unconstrained vector in \((-\infty, \infty)\)

Return Value: vector A vector of coefficients for a stationary AR or inverible MA process.

Constrain vector of coefficients to the stationary and intertible region for AR or MA functions.

See R. H. Jones (1980), M. C. Jones (1987), Monahan (1984), Ansley and Kohn (1986), and the functions tools.constrain_stationary_univariate and tools.unconstraine_stationary_univariate in statsmodels.tsa.statespace.

vector constrain_stationary(vector x) {
  vector[num_elements(x)] r;
  int n;
  n = num_elements(x);
  // transform (-Inf, Inf) to (-1, 1)
  for (i in 1:n) {
    r[i] = x[i] / (sqrt(1.0 + pow(x[i], 2)));
  }
  // Transform PACF to ACF
  return pacf_to_acf(r);
}

4.9.3 acf_to_pacf

Parameters:

• vector x Coeffcients of an autocorrelation function.

Return Value: vector A vector of coefficients of the corresponding partial autocorrelation function.

Convert coefficients of an autocorrelation function to partial autocorrelations.

vector acf_to_pacf(vector x) {
  matrix[num_elements(x), num_elements(x)] y;
  vector[num_elements(x)] r;
  int n;
  n = num_elements(x);
  y = rep_matrix(0.0, n, n);
  y[n] = -x ';
  for (j in 0:(n - 1)) {
    int k;
    k = n - j;
    for (i in 1:(k - 1)) {
      y[k - 1, i] = (y[k, i] - y[k, k] * y[k, k - i]) / (1 - pow(y[k, k], 2));
    }
  }
  r = diagonal(y);
  return r;
}

4.9.4 unconstrain_stationary

Parameters:

• vector x Coeffcients of an autocorrelation function.

Return Value: vector Coefficients of the corresponding partial autocorrelation function.

Transform from stationary and invertible space to \((-\infty, \infty)\).

vector unconstrain_stationary(vector x) {
  matrix[num_elements(x), num_elements(x)] y;
  vector[num_elements(x)] r;
  vector[num_elements(x)] z;
  int n;
  n = num_elements(x);
  // Transform ACF to PACF
  r = acf_to_pacf(x);
  // Transform (-1, 1) to (-Inf, Inf)
  for (i in 1:n) {
    z[i] = r[i] / (sqrt(1.0 - pow(r[i], 2)));
  }
  return z;
}

4.9.5 kronecker_prod

Parameters:

• matrix A An \(m \times n\) matrix
• matrix B A \(p \times q\) matrix

Return Value: matrix An \(mp \times nq\) matrix.

Kronecker product

The Kronecker product of a \(A\) and \(B\) is \[ A \otimes B = \begin{bmatrix} a_{11} B \cdots a_{1n} B \\ \vdots & \ddots & vdots \\ a_{m1} B & \cdots & a_{mn} B \end{bmatrix} . \]

matrix kronecker_prod(matrix A, matrix B) {
  matrix[rows(A) * rows(B), cols(A) * cols(B)] C;
  int m;
  int n;
  int p;
  int q;
  m = rows(A);
  n = cols(A);
  p = rows(B);
  q = cols(B);
  for (i in 1:m) {
    for (j in 1:n) {
      int row_start;
      int row_end;
      int col_start;
      int col_end;
      row_start = (i - 1) * p + 1;
      row_end = (i - 1) * p + p;
      col_start = (j - 1) * q + 1;
      col_end = (j - 1) * q + 1;
      C[row_start:row_end, col_start:col_end] = A[i, j] * B;
    }
  }
  return C;
}

4.9.6 arima_stationary_cov

Parameters:

• matrix T The \(m \times m\) transition matrix
• matrix R The \(m \times q\) system disturbance selection matrix

Return Value: matrix An \(m \times m\) matrix with the stationary covariance matrix.

Find the covariance of the stationary distribution of an ARMA model

The initial conditions are \(\alpha_1 \sim N(0, \sigma^2 Q_0)\), where \(Q_0\) is the solution to \[ (T \otimes T) \VEC(Q_0) = \VEC(R R') \] where \(\VEC(Q_0)\) and \(\VEC(R R')\) are the stacked columns of \(Q_0\) and \(R R'\)

See Durbin and Koopman (2012), Sec 5.6.2.

matrix arima_stationary_cov(matrix T, matrix R) {
  matrix[rows(T), cols(T)] Q0;
  matrix[rows(T) * rows(T), rows(T) * rows(T)] TT;
  vector[rows(T) * rows(T)] RR;
  int m;
  int m2;
  m = rows(T);
  m2 = m * m;
  RR = to_vector(tcrossprod(R));
  TT = kronecker_prod(T, T);
  Q0 = to_matrix_colwise((diag_matrix(rep_vector(1.0, m2)) - TT) \ RR, m, m);
  return Q0;
}

References

Durbin, J., and S.J. Koopman. 2012. Time Series Analysis by State Space Methods: Second Edition. Oxford Statistical Science Series. OUP Oxford. http://books.google.com/books?id=fOq39Zh0olQC.

Durbin, J., and S. J. Koopman. 2002. “A Simple and Efficient Simulation Smoother for State Space Time Series Analysis.” Biometrika 89 (3). Biometrika Trust: 603–15. http://www.jstor.org/stable/4140605.

Jones, Richard H. 1980. “Maximum Likelihood Fitting of ARMA Models to Time Series with Missing Observations.” Technometrics 22 (3). [Taylor & Francis, Ltd., American Statistical Association, American Society for Quality]: 389–95. http://www.jstor.org/stable/1268324.

Jones, M. C. 1987. “Randomly Choosing Parameters from the Stationarity and Invertibility Region of Autoregressive-Moving Average Models.” Journal of the Royal Statistical Society. Series C (Applied Statistics) 36 (2). [Wiley, Royal Statistical Society]: 134–38. http://www.jstor.org/stable/2347544.

Monahan, John F. 1984. “A Note on Enforcing Stationarity in Autoregressive-Moving Average Models.” Biometrika 71 (2): 403–4. doi:10.1093/biomet/71.2.403.

Ansley, Craig F., and Robert Kohn. 1986. “A Note on Reparameterizing a Vector Autoregressive Moving Average Model to Enforce Stationarity.” Journal of Statistical Computation and Simulation 24 (2): 99–106. doi:10.1080/00949658608810893.