statespace_simsmooth#

Simulation from the model and the simulation smoothers (DK §4.9, §5.5).

The random input is standard normal variates, so draws are reproducible: u_init (m), u_eps (p, n) and u_eta (r, n). They are scaled by lower triangular square roots of \(P_*\), \(H_t\) and \(Q_t\) (the Cholesky factor when the matrix is positive definite), as in statsmodels, so the same variates give the same draws.

simsmooth_result_t#

type :: simsmooth_result_t
  real(dp), allocatable :: state(:, :)  ! (m, n) draw of alpha_t | Y_n
  real(dp), allocatable :: eps(:, :)    ! (p, n) draw of eps_t | Y_n
  real(dp), allocatable :: eta(:, :)    ! (r, n) draw of eta_t | Y_n
end type

simulate#

subroutine simulate(rep, u_init, u_eps, u_eta, y, alpha, eps, eta, info)
  real(dp), intent(in) :: u_init(:), u_eps(:, :), u_eta(:, :)
  real(dp), intent(out) :: y(:, :), alpha(:, :), eps(:, :), eta(:, :)

Simulate from the model:

\[\alpha_1 = a_1 + S(P_*) u_{init}, \quad y_t = d_t + Z_t \alpha_t + S(H_t) u_{\varepsilon,t}, \quad \alpha_{t+1} = c_t + T_t \alpha_t + R_t S(Q_t) u_{\eta,t},\]

with S(A) a lower triangular square root of A. Missing values in rep%y are ignored. Diffuse states start at their mean.

simulation_smoother#

subroutine simulation_smoother(rep, u_init, u_eps, u_eta, sim, info)
  type(simsmooth_result_t), intent(out) :: sim

Draw \((\alpha, \varepsilon, \eta)\) given the data by mean corrections (Durbin and Koopman 2002; DK §4.9.2):

  1. simulate \((y^+, \alpha^+, \varepsilon^+, \eta^+)\) from the model;

  2. smooth \(y - y^+\) in the model with zero intercepts and \(a_1 = 0\), which gives \(\hat\alpha(y) - \hat\alpha(y^+)\) because smoothing is affine;

  3. add the result to the simulated values.

Diffuse directions of \(\alpha_1\) get no draw: the exact diffuse smoother reproduces any shift in them, so they cancel in step 3 (DK §5.5). For a missing element the drawn \(\varepsilon\) is conditional on the observed elements (statsmodels draws it unconditionally).

djs_measurement_disturbances#

subroutine djs_measurement_disturbances(rep, fres, u_eps, eps, info)
  type(filter_result_t), intent(in) :: fres
  real(dp), intent(in) :: u_eps(:, :)
  real(dp), intent(out) :: eps(:, :)

The de Jong-Shephard simulation smoother for \(\varepsilon\) (DK §4.9.3, eqs. 4.83-4.88): going back from t = n, draw \(\varepsilon_t\) from its distribution given the data and the later draws. Not defined with diffuse states.

djs_state_disturbances#

subroutine djs_state_disturbances(rep, fres, u_eta, u_init, eta, alpha, info)
  real(dp), intent(in) :: u_eta(:, :), u_init(:)
  real(dp), intent(out) :: eta(:, :), alpha(:, :)

The same for \(\eta\) (DK eqs. 4.89-4.91), then the states forward (DK eq. 4.92). DK start the states from the mean \(a_1 + P_1 \tilde r_0\) of \(\alpha_1\) given the data and the \(\eta\) draws; this routine draws \(\alpha_1\) from that distribution, with variance \(P_1 - P_1 \tilde N_0 P_1\), using u_init. Without that draw the state draws are too concentrated. Not defined with diffuse states.

psd_sqrt#

function psd_sqrt(A) result(S)

Lower triangular S with \(S S' = A\) for symmetric positive semi-definite A: \(S = L D^{1/2}\) from \(A = L D L'\); the Cholesky factor when A is positive definite.

draw_standard_normal#

subroutine draw_standard_normal(x)
  real(dp), intent(out) :: x(:)

Fill x with independent N(0, 1) draws (Box-Muller on random_number).