WIP: Exogenous inputs + EM estimation

[GordStephen]

I wouldn't want to pull this yet but I though I'd let people know what I've been up to - implementing deterministic inputs to the evolution and observation equations, with appropriate EM estimation.

I added a new type called StateSpaceModelMask that mirrors a SSM and lets you define which of the model parameters you would want to estimate. Unfixed parameter estimation (estimate every variable possible) should be working, as should estimation with control input parameters fixed (e.g. to zeros). Initial results are encouraging relative to the current estimator - using the current fit test (on a longer generated series):

Current (Nelder-Mead / Optim.jl)

julia> @time fit(y, build, randn(9))
 36.789447 seconds (323.65 M allocations: 14.849 GB, 6.64% gc time)
([1.000598293008593,-10.99513419805376,-0.20039712192006387,-0.4011806057563814,0.0997647355757492,7.218337983064723,8.552265107672367,5.801225460724416,-12.392217275096918],StateSpaceModel{Float64}, 1-D process x 3-D observations
Process evolution matrix F:
[1.000598293008593]

Process error covariance V:
[1.678316599450008e-5]

Observation matrix G:
[-0.20039712192006387
 -0.4011806057563814
 0.0997647355757492]

Obseration error covariance W:
[1364.2198190288902 0.0 0.0
 0.0 5178.4709468500205 0.0
 0.0 0.0 330.70457716329184],KalmanSmoothed{Float64}
1000 observations, 1-D process x 3-D observations
Negative log-likelihood: 13651.098608054745
)

EM:

julia> @time fit(y, build(randn(9)))
  1.793928 seconds (6.30 M allocations: 731.435 MB, 5.24% gc time)
(StateSpaceModel{Float64}, 1-D process x 3-D observations
Process evolution matrix F:
[0.9994321799454007]

Process error covariance V:
[0.7466870359836347]

Observation matrix G:
[1.830905476396109
 3.6618367191926926
 -0.91132011739738]

Obseration error covariance W:
[7.658804608602309 0.0 0.0
 0.0 2.18631683886998 0.0
 0.0 0.0 3.8273984210163285],KalmanSmoothed{Float64}
1000 observations, 1-D process x 3-D observations
Negative log-likelihood: 7260.434009339921

And for reference, the model that originally generated the series is

julia> mod1
StateSpaceModel{Float64}, 1-D process x 3-D observations
Process evolution matrix F:
1x1 Array{Float64,2}:
 1.0

Process error covariance V:
1x1 Array{Float64,2}:
 2.0

Observation matrix G:
3x1 Array{Float64,2}:
  1.0
  2.0
 -0.5

Obseration error covariance W:
3x3 Array{Float64,2}:
 8.0  0.0  0.0
 0.0  2.5  0.0
 0.0  0.0  4.0

So, significantly faster and better results. Yay :)

I also just found this, which lays things out explicitly and would would have saved me a lot of time on the weekend! https://cran.r-project.org/web/packages/MARSS/vignettes/EMDerivation.pdf Next step is implementing the constrained estimation techniques it describes.

[milktrader]

Sweet!

[GordStephen]

More progress. Added new types to support the specification of linear-constrained parameter matrices (e.g. [a b; b 2a]) and implemented EM estimation for those constraints. Two bigger outstanding issues are 1) support for fitting parameters to deterministic processes (chapters 7-10 in the link above - eg zeros in the noise covariance matrix diagonal - right now the fitter needs to invert those matrices which is an issue) and 2) I took out time-dependent transition/observation matrix support and haven't added it back. I'll get the Julia v0.3 tests passing at some point as well.

[milktrader]

Do you want this PR to work on Julia 0.3? I'm restricting the package to work only on Julia 0.3 for the time being and then once we have that landscape cleaned up we can move development to using Julia 0.4 only.

UPDATE: tag 0.0.3 will be the first step to Julia 0.3 isolation. Once we begin on Julia 0.4, a minor tag will be used (i.e., 0.1.0)

[GordStephen]

No, I'm fine for this to be 0.4-only.

[milktrader]

Ok, I'm soon tagging 0.0.3, which will include support for both Julia 0.3 and Julia 0.4. We've gotten into supporting both versions of Julia too much to isolate at this point. The 0.1.0 tag will support only Julia 0.4 though, but we're still on release candidates.

[GordStephen]

Changes here were either introduced in #52 or will be in #57.