The statespacemodels.jl from lampspuc

There is a missing minus sign in the seasonality equation.

possible test failure in upcoming Julia version 1.5

A PkgEval run for a PR which changes the generated numbers for randn! indicates that the tests of this package might fail in Julia 1.5 (and on Julia current master). Apologies if this is a false positive.
cf. https://github.com/JuliaCI/NanosoldierReports/blob/7de24e455342298cbef56826b5827f0d7640d2c1/pkgeval/by_hash/b89e35c_vs_098ef24/logs/StateSpaceModels/1.5.0-DEV-71a4a114c2.log

Add time limit option to optimization

Optim.Options has time_limit optition. It would be nice to add this to LBFGS struct.

Make an example where there are no hyperparameters

We should write an example for the usage of the online filter where there is no need for estimation.

Change structures to Array{Type, 3}

This code

        y = collect(1:15.0)
        unimodel = lineartrendmodel(y)

        @test isa(unimodel, StateSpaceModels.StateSpaceModel)
        @test unimodel.mode == "time-invariant"

        ss = statespace(unimodel)

        @test isa(ss, StateSpaceModels.StateSpace)
        ss.state.alpha

returns this as output, it is a super annoying structure to deal with

UnivariateKalmanFilter failing for log airline

The function forecast should return a forecast not an array of distributions

Improve documentation / docstrings

Change Kronecker product order when building system matrices

Add examples

Section 8.6 from Koopman book has an example of Diebold-Li
Examples from PyFlux to compare results: Nile River, DAR, US log-GDP
Example with NaN to estimate a parameter in Z or T

Adjust steady state check for time-variant case

Add diagnostics struct to be returned after model is estimated

Implement exact initialization filter

Implement EM algorithm for MLE

Decide how to model Z[t] for t > n for forecasting

Add `verbose` key argument

I think there should be a verbose key argument to make progress logs optional.

Add Plot recipes for forecasting and components

Here we could create another abstraction such as an UnobservedComponentsModel where we should give a name to each component

Handle observation missing values

Add AIC and BIC as diagnostics for estimation

Update airline example in docs to the version in examples.jl

Question about long structural cycles

Hi,

Thanks for the really cool module.

I have data that has long "seasonal" cycles, think of data being sampled each minute and each day shows cyclical changes. It seems that the larger I make the s parameter of the structural model the longer and longer it takes to fit.

For example, it can take a really long time for 1000 data-points and 288 specified for s. Is there some trick I can use to make things go faster?

Thanks,

Rusty

Add HypothesisTest to run Diagnostics

Should we consider using HypothesisTests instead of creating our own version? This could be at least on the tests to assure we are doing the right thing.

Documentation of estimate function is out of data

The function should be statespace(model; filter_type = KalmanFilter, opt_method = LBFGS(), verbose = 1)

To do list

Implement diagnostics

Diagnostics for the residuals:

Jarque-Bera test
Ljung-Box test
Homoscedasticity test

Conflicting method definition

Following method definition conflict occur.

Method definition vec(Number) in module FiniteDiff at /home/appleparan/.julia/packages/FiniteDiff/IPtYL/src/jacobians.jl:128 overwritten in module DiffEqDiffTools at /home/appleparan/.julia/packages/DiffEqDiffTools/3mm8U/src/jacobians.jl:114.
  ** incremental compilation may be fatally broken for this module **

It seems Optim package version limit in Project.toml makes this warning. Is there a reason keeping version limit in Project.toml?

PS: here is a dependencies of problematic packages.

Optim -> NLSolversBase -> FiniteDiff
Optim -> DiffEqDiffTools

Eliminate StateSpaceCovariance struct

This is one of the first steps to make a more general estimation

Function to update just the Kalman filter

Maybe it is good to have a function to just update the Kalman filter. My sketch:

function update_kalman_filter(model::StateSpaceModel{Typ}, filter::FilterOutput, model_test::StateSpaceModel{Typ}; tol::Typ = Typ(1e-5)) where Typ
    time_invariant = model.mode == "time-invariant"

    # Predictive state and its covariance
    a = cat(filter.a, Matrix{Typ}(undef, model_test.dim.n, model.dim.m), dims=1)
    P = cat(filter.P, Array{Typ, 3}(undef, model.dim.m, model.dim.m, model_test.dim.n+1), dims=3)
    att = cat(filter.att, Matrix{Typ}(undef, model_test.dim.n, model.dim.m), dims=1)
    Ptt = cat(filter.Ptt, Array{Typ, 3}(undef, model.dim.m, model.dim.m, model_test.dim.n), dims=3)

    # Innovation and its covariance
    v = cat(filter.v, Matrix{Typ}(undef, model_test.dim.n, model.dim.p), dims=1)
    F = cat(filter.F, Array{Typ, 3}(undef, model.dim.p, model.dim.p, model_test.dim.n) , dims=3)

    # Kalman gain
    K = Array{Typ, 3}(undef, model.dim.m, model.dim.p, model.dim.n + model_test.dim.n)
    
    # Auxiliary structures
    ZP = Array{Typ, 2}(undef, model.dim.p, model.dim.m)
    P_Ztransp_invF = Array{Typ, 2}(undef, model.dim.m, model.dim.p)
    RQR = Array{Typ, 2}(undef, model.dim.m, model.dim.m)

    # Steady state initialization
    steadystate = filter.steadystate
    tsteady     = filter.tsteady

    # Initial state: big Kappa initialization
    StateSpaceModels.fill_a1!(a)
    StateSpaceModels.fill_P1!(P; bigkappa = Typ(1e6))
    y = vcat(model.y, model_test.y)
    d = vcat(model.d, model_test.d)
    Z = cat(model.Z, model_test.Z,dims=3)
    c = vcat(model.c, model_test.c)
    
    StateSpaceModels.update_ZP!(ZP, Z, P, filter.tsteady) # ZP = Z[:, :, t] * P[:, :, t]
    StateSpaceModels.update_F!(F, ZP, Z, model.H, filter.tsteady) # F_t = Z_t * P_t * Z_t + H
    StateSpaceModels.update_P_Ztransp_Finv!(P_Ztransp_invF, ZP, F, filter.tsteady) # P_Ztransp_invF   = ZP' * invF(F, t)
    StateSpaceModels.update_K!(K, P_Ztransp_invF, model.T, filter.tsteady) # K_t = T * P_t * Z_t * F^-1_t
    
    for t_ = 1:model_test.dim.n
        t = model.dim.n+t_-1
        if t in model_test.missing_observations
            steadystate  = false
            v[t, :]      = fill(NaN, model.dim.p)
            F[:, :, t]   = fill(NaN, (model.dim.p, model.dim.p))
            StateSpaceModels.update_att!(att, a, t) # attt = a_t
            StateSpaceModels.update_Ptt!(Ptt, P, t) # Pttt = P_t
            StateSpaceModels.update_a!(a, att, model.T, c, t) # a_t+1 = T * attt + c_t
            StateSpaceModels.update_P!(P, model.T, Ptt, RQR, t) # P_t+1 = T * Pttt * T' + RQR'
        elseif steadystate
            StateSpaceModels.update_v!(v, y, Z, d, a, t) # v_t = y_t - Z_t * a_t - d_t
            StateSpaceModels.repeat_matrix_t_plus_1!(F, t-1) # F[:, :, t]   = F[:, :, t-1]
            StateSpaceModels.repeat_matrix_t_plus_1!(K, t-1) # K[:, :, t]   = K[:, :, t-1]
            StateSpaceModels.update_att!(att, a, P_Ztransp_invF, v, t) # attt = a_t + P_t * Z_t * F^-1_t * v_t
            StateSpaceModels.repeat_matrix_t_plus_1!(Ptt, t-1) # Pttt = Pttt-1
            StateSpaceModels.update_a!(a, att, model.T, c, t) # a_t+1 = T * attt + c_t
            StateSpaceModels.repeat_matrix_t_plus_1!(P, t) # P__+1 = P_t
        else
            StateSpaceModels.update_v!(v, y, Z, d, a, t) # v_t = y_t - Z_t * a_t - d_t
            StateSpaceModels.update_ZP!(ZP, Z, P, t) # ZP = Z[:, :, t] * P[:, :, t]
            StateSpaceModels.update_F!(F, ZP, Z, model.H, t) # F_t = Z_t * P_t * Z_t + H
            StateSpaceModels.update_P_Ztransp_Finv!(P_Ztransp_invF, ZP, F, t) # P_Ztransp_invF   = ZP' * invF(F, t)
            StateSpaceModels.update_K!(K, P_Ztransp_invF, model.T, t) # K_t = T * P_t * Z_t * F^-1_t
            StateSpaceModels.update_att!(att, a, P_Ztransp_invF, v, t) # attt = a_t + P_t * Z_t * F^-1_t * v_t
            StateSpaceModels.update_Ptt!(Ptt, P, P_Ztransp_invF, ZP, t) # Pttt = P_t - P_t * Z_t' * F^-1_t * Z_t * P_t
            StateSpaceModels.update_a!(a, att, model.T, c, t) # a_t+1 = T * attt + c_t
            StateSpaceModels.update_P!(P, model.T, Ptt, RQR, t) # P_t+1 = T * Pttt * T' + RQR'
            if time_invariant && StateSpaceModels.check_steady_state(P, t, tol)
                steadystate = true
                tsteady     = t
            end
        end
    end
    # Return the auxiliary filter structre
    return KalmanFilter(a[model.dim.n+1:end,:], att[model.dim.n+1:end,:], v, P, Ptt, F, steadystate, tsteady, K)
end

Implement standard Kalman filter

Shift between forecast and simulation scenarions in 0.4 branch

From the "Basic Structural Model" test in the refactor-to-0.4 branch:

julia> forec.expected_value
10-element Array{Array{Float64,1},1}:
 [6.083167302258863]
 [6.1946408652357725]
 [6.2159346905082975]
 [6.224800298528862]
 [6.342665144161216]
 [6.478336552034039]
 [6.4752275105634185]
 [6.305246285595047]
 [6.2049769753007]
 [6.068300933194807]

julia> expected_value_of_scenarios(scenarios)
10-element Array{Array{Float64,1},1}:
 [6.125259198552738]
 [6.082992264095608]
 [6.19469004566561]
 [6.215805308795272]
 [6.224534124065003]
 [6.342838593417967]
 [6.478180596529146]
 [6.475235575886917]
 [6.305295124889548]
 [6.205057985402772]

There is a clear shift of one time period between both.

Note that this should be caught in the tests, but they're currently using atol = 1.0 which is super large. We should use rtol = 1e-3 or even lower for this.

Update docs

Docs are currently out of date with the new change allowing any float type for the models and filters

Forecast and Simulate are missing d and c matrices

Add ARIMA constructor

Potential assertion error

These functions are missing @assertion macros to avoid matricesandcto have different lenghts fromy`

    function StateSpaceModel(y::VecOrMat{Typ}, Z::Array{Typ, 3}, T::Matrix{Typ}, R::Matrix{Typ}, 
                             d::Matrix{Typ}, c::Matrix{Typ}, H::Matrix{Typ}, Q::Matrix{Typ}) where Typ <: Real

        # Convert y to a Matrix
        y = ensure_is_matrix(y)
        # Build StateSpaceDimensions
        dim = build_ss_dim(y, Z, T, R)
        return new{Typ}(y, Z, T, R, d, c, H, Q, dim, find_missing_observations(y), "time-variant")
    end

    function StateSpaceModel(y::VecOrMat{Typ}, Z::Matrix{Typ}, T::Matrix{Typ}, R::Matrix{Typ}, 
                             d::Matrix{Typ}, c::Matrix{Typ}, H::Matrix{Typ}, Q::Matrix{Typ}) where Typ <: Real

        # Convert y to a Matrix
        y = ensure_is_matrix(y)
        # Build StateSpaceDimensions
        dim = build_ss_dim(y, Z, T, R)
        Zvar = Array{Typ, 3}(undef, dim.p, dim.m, dim.n)
        for t in 1:dim.n, i in axes(Z, 1), j in axes(Z, 2)
            Zvar[i, j, t] = Z[i, j]
        end

        return new{Typ}(y, Zvar, T, R, d, c, H, Q, dim, find_missing_observations(y), "time-invariant")
    end

lampspuc / statespacemodels.jl Goto Github PK

statespacemodels.jl's People

Contributors

Stargazers

Watchers

Forkers

statespacemodels.jl's Issues

Recommend Projects

Recommend Topics

Recommend Org