statespace_mle

Maximum likelihood estimation (DK §7.3) with L-BFGS-B (Zhu, Byrd, Lu and Nocedal 1997; Morales and Nocedal 2011).

The optimizer minimizes -loglike / n over the unconstrained parameters, as statsmodels does, with statsmodels’ default tolerances. The gradient is the analytic score (DK §7.3.3; statespace_score) for the parameters it covers and central differences for the others (parameters in Z or T, as DK recommend). Standard errors come from a numerical Hessian in the unconstrained parameters, mapped to the constrained ones by the delta method (DK §7.3.6). Points where the likelihood cannot be evaluated get a large objective value, so the line search backs off.

See Estimation for the choices behind the gradient and the standard errors.

Example

type(fit_result_t) :: res
type(fit_options_t) :: opts

opts%factr = 10.0_dp        ! tight convergence
opts%pgtol = 1.0e-9_dp
call fit(model, res, options=opts, info=info)
print *, res%params, res%bse, res%llf

fit_options_t

type :: fit_options_t
  integer :: maxiter = 500
  integer :: m = 10                  ! L-BFGS corrections
  real(dp) :: factr = 1.0e7_dp       ! relative reduction tolerance / machine epsilon
  real(dp) :: pgtol = 1.0e-5_dp      ! projected gradient tolerance
  logical :: compute_cov = .true.
  integer :: iprint = -1             ! L-BFGS-B output; negative is silent
  integer :: gradient = GRADIENT_AUTO
end type

gradient is GRADIENT_AUTO (0; analytic where it applies), GRADIENT_NUMERICAL (1) or GRADIENT_ANALYTIC (2; fail with SS_ERR_UNSUPPORTED if no parameter is covered).

fit_result_t

type :: fit_result_t
  real(dp), allocatable :: params(:)          ! constrained estimates
  real(dp), allocatable :: cov_params(:, :)
  real(dp), allocatable :: bse(:)
  real(dp) :: llf, scale, aic, bic
  integer :: niter, nfev
  logical :: converged, analytic_gradient
  character(len=60) :: message
end type

aic and bic count the diffuse states and a concentrated scale as parameters and use n minus the burn-in, as statsmodels does; DK §7.4 also divide by n. analytic_gradient is true if the score covered at least one parameter.

fit

subroutine fit(model, res, start_params, options, info)
  class(ssm_model_t), intent(inout) :: model
  type(fit_result_t), intent(out) :: res
  real(dp), intent(in), optional :: start_params(:)
  type(fit_options_t), intent(in), optional :: options
  integer, intent(out) :: info

Estimate the parameters; the model is left at the estimates. If the estimates exist but the Hessian is not negative definite, info is SS_ERR_NOT_PD and cov_params and bse are not allocated.

fit_many

subroutine fit_many(models, res, info, options)
  class(ssm_model_t), intent(inout) :: models(:)
  type(fit_result_t), intent(out) :: res(:)
  integer, intent(out) :: info(:)

Fit independent models in parallel with OpenMP when the library is built with -fopenmp, in turn otherwise. Each model is touched by one thread only; L-BFGS-B keeps its state in the caller’s arrays, so it is reentrant. For large models set OPENBLAS_NUM_THREADS=1 to avoid nested threading.

numerical_hessian

function numerical_hessian(model, params, info, transformed) result(hess)

Central-difference Hessian of the log likelihood at params, constrained unless transformed = .false..

estimation_bias

subroutine estimation_bias(model, fres, ndraw, bias_alpha, info, bias_V, antithetic, failed)
  type(fit_result_t), intent(in) :: fres
  integer, intent(in) :: ndraw
  real(dp), intent(out) :: bias_alpha(:, :)                ! (m, n)
  real(dp), intent(out), optional :: bias_V(:, :, :)       ! (m, m, n)
  logical, intent(in), optional :: antithetic              ! default .true.
  integer, intent(out), optional :: failed

The bias of the smoothed state from estimating the parameters (DK §7.3.7, eq. 7.20). Treating \(\hat\psi\) as the true value, draw \(\psi^{(i)} \sim N(\hat\psi, \Omega)\) on the unconstrained scale, with \(\Omega\) from fres%cov_params by the delta method, and estimate

\[B = \frac1N \sum_i \hat\alpha(\psi^{(i)}) - \hat\alpha(\hat\psi).\]

With antithetic, each draw is paired with its reflection \(2\hat\psi - \psi^{(i)}\), and ndraw must be even. Draws where the model cannot be evaluated are skipped with their pair and counted in failed. Uses random_number; seed it for reproducible results. The model is left at \(\hat\psi\).