Quickstart#

The local level model for the flow of the Nile (DK ch. 2),

\[y_t = \mu_t + \varepsilon_t, \quad \varepsilon_t \sim N(0, \sigma^2_\varepsilon), \qquad \mu_{t+1} = \mu_t + \eta_t, \quad \eta_t \sim N(0, \sigma^2_\eta),\]

estimated by maximum likelihood with an exact diffuse initial level. DK report \(\hat\sigma^2_\varepsilon = 15099\) and \(\hat\sigma^2_\eta = 1469.1\); the likelihood is flat near its maximum, so optimizers agree to about four digits.

Python#

A built-in model:

>>> import numpy as np
>>> import ssfortran as ss
>>> y = np.loadtxt("data/nile.csv", delimiter=",", skiprows=1)[:, 1]
>>> mod = ss.StructuralModel(y, [ss.Irregular(), ss.Level()])
>>> res = mod.fit(factr=10, pgtol=1e-9)
>>> np.round(res.params)
array([15099.,  1469.])
>>> round(res.llf, 4)
-633.4646

The smoothed level and its variance:

>>> sm = res.smooth()
>>> level = sm.smoothed_state[0]
>>> np.round(level[:3])
array([1112., 1111., 1105.])

The same model declared from its matrices, with the two variances as parameters (Defining models):

>>> mod = ss.MappedModel(y, k_states=1, k_params=2,
...                      param_names=["sigma2.irregular", "sigma2.level"],
...                      start_params=[np.var(y) / 2, np.var(y) / 2])
>>> mod["design"] = mod["transition"] = mod["selection"] = [[1.0]]
>>> mod.initialize_diffuse()
>>> _ = mod.map(0, "obs_cov", 0, 0).map(1, "state_cov", 0, 0)
>>> _ = mod.constrain([0, 1], "positive")
>>> res = mod.fit(factr=10, pgtol=1e-9)
>>> np.round(res.params)
array([15099.,  1469.])

Every model builds its Representation, the system at given parameter values. mod.representation(res.params) returns it, to run the other algorithms of DK Part I (Filtering and smoothing).

Fortran#

The same model defined by extending ssm_model_t, from example/nile_mle.f90:

program nile_mle
  use statespace
  use nile_local_level_model, only: local_level_t, local_level
  implicit none

  type(local_level_t) :: model
  type(fit_result_t) :: res
  type(fit_options_t) :: opts
  character(len=32), allocatable :: names(:)
  real(dp), allocatable :: y(:, :)
  integer :: info, i

  y = read_nile("data/nile.csv")
  model = local_level(y)

  opts%factr = 10.0_dp       ! tight convergence; the likelihood is flat
  opts%pgtol = 1.0e-9_dp
  call fit(model, res, options=opts, info=info)
  if (info /= SS_OK) error stop "fit failed"

  names = model%param_names()
  print '(a)', trim(res%message)
  print '(a, i0, a, i0)', "iterations: ", res%niter, "  likelihood evaluations: ", &
      res%nfev
  print '(a, f14.6)', "log likelihood: ", res%llf
  print '(a, f10.4, a, f10.4)', "AIC: ", res%aic, "  BIC: ", res%bic
  print '(a20, 2a14)', "", "estimate", "std err"
  do i = 1, model%k_params
    print '(a20, 2f14.4)', trim(names(i)), res%params(i), res%bse(i)
  end do

contains

  function read_nile(path) result(y)
    character(len=*), intent(in) :: path
    real(dp), allocatable :: y(:, :)
    real(dp) :: year, vol
    integer :: unit, ios, n

    open (newunit=unit, file=path, status="old", action="read")
    read (unit, *)
    n = 0
    do
      read (unit, *, iostat=ios) year, vol
      if (ios /= 0) exit
      n = n + 1
    end do
    rewind (unit)
    read (unit, *)
    allocate (y(1, n))
    do n = 1, size(y, 2)
      read (unit, *) year, y(1, n)
    end do
    close (unit)
  end function read_nile
end program nile_mle

Output:

$ fpm run --profile release --example nile_mle
CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL
iterations: 10  likelihood evaluations: 12
log likelihood:    -633.464564
AIC:  1272.9291  BIC:  1280.7446
                          estimate       std err
    sigma2.irregular    15098.5183     3145.5481
        sigma2.level     1469.1764     1280.3752

Next#

  • User guide explains the model, the algorithms and the choices behind them.

  • Examples reproduces the illustrations of DK ch. 8.

  • Reference lists every class, routine and argument.