first steps with filtering package

Project.toml

@@ -6,9 +6,11 @@ version = "0.1.29"
 [deps]
 AbstractDifferentiation = "c29ec348-61ec-40c8-8164-b8c60e9d9f3d"
 AxisKeys = "94b1ba4f-4ee9-5380-92f1-94cde586c3c5"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 BlockTriangularForm = "adeb47b7-70bf-415a-bb24-c358563e873a"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 DynarePreprocessor_jll = "23afba7c-24e5-5ee2-bc2c-b42e07f0492a"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
@@ -18,6 +20,7 @@ Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearOperators = "5c8ed15e-5a4c-59e4-a42b-c7e8811fb125"
+LowLevelParticleFilters = "d9d29d28-c116-5dba-9239-57a5fe23875b"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 REPL = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
@@ -29,6 +32,7 @@ RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 SpeedMapping = "f1835b91-879b-4a3f-a438-e4baacf14412"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Subscripts = "2b7f82d5-8785-4f63-971e-f18ddbeb808e"
 SymPyPythonCall = "bc8888f7-b21e-4b7c-a06a-5d9c9496438c"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"

test/figure_out_filtering.jl

@@ -0,0 +1,103 @@
+using MacroModelling, LowLevelParticleFilters, Distributions
+import LinearAlgebra as ℒ
+
+include("models/FS2000.jl")
+
+simulations = get_simulation(m)
+
+# (x,u,p,t) = (state, input, parameters, time)
+# state_update₁ = function(state::Vector{Float64}, shock::Vector{Float64}) sol_mat * [state[𝓂.timings.past_not_future_and_mixed_idx]; shock] end
+
+# dynamics(state, input, parameters, time) = m.solution.perturbation.first_order.solution_matrix * [state[m.timings.past_not_future_and_mixed_idx]; input]
+dynamics(state, input, parameters, time) = m.solution.perturbation.first_order.state_update(state,input)
+measurement(state, input, parameters, time) = m.solution.perturbation.first_order.state_update(state,input)#ℒ.I * state
+
+covariance = get_covariance(m) |> collect
+# covariance = sparse(covariance)
+# droptol!(covariance,eps())
+# covariance = sparse(Symmetric(covariance))
+# isposdef(covariance)
+ukf = UnscentedKalmanFilter(dynamics, 
+                            measurement, 
+                            eye(m.timings.nVars), 
+                            eye(m.timings.nVars), 
+                            MvNormal(diag(covariance)), 
+                            nu = m.timings.nExo, 
+                            ny = m.timings.nVars)
+
+
+x,u,y = LowLevelParticleFilters.simulate(ukf,10,MvNormal(ones(m.timings.nExo))) # Simuate trajectory using the model in the filter
+
+
+using LinearAlgebra
+covariance
+isposdef((covariance+covariance')/2)
+factorize(Symmetric(covariance))
+
+
+
+
+using Statistics, LinearAlgebra, StaticArrays
+m = randn(3)
+S = randn(3,3)
+S = S'S
+xs = LowLevelParticleFilters.sigmapoints(m, S)
+X = reduce(hcat, xs)
+@test vec(mean(X, dims=2)) ≈ m
+@test Statistics.cov(X, dims=2) ≈ S
+
+m = [1,2]
+S = [3. 1; 1 4]
+xs = LowLevelParticleFilters.sigmapoints(m, S)
+X = reduce(hcat, xs)
+# @test vec(mean(X, dims=2)) ≈ m
+# @test cov(X, dims=2) ≈ S
+
+
+
+
+eye(n) = Matrix{Float64}(I,n,n)
+mvnormal(d::Int, σ::Real) = MvNormal(LinearAlgebra.Diagonal(fill(float(σ) ^ 2, d)))
+mvnormal(μ::AbstractVector{<:Real}, σ::Real) = MvNormal(μ, float(σ) ^ 2 * I)
+
+nx = 2 # Dinemsion of state
+nu = 2 # Dinemsion of input
+ny = 2 # Dinemsion of measurements
+
+d0 = mvnormal(randn(nx),2.0)   # Initial state Distribution
+du = mvnormal(2,1) # Control input distribution
+
+# Define random linenar state-space system
+Tr = randn(nx,nx)
+A = SA[0.99 0.1; 0 0.2]
+B = @SMatrix randn(nx,nu)
+C = SMatrix{ny,ny}(eye(ny))
+# C = SMatrix{p,n}([1 1])
+
+dynamics(x,u,p,t) = A*x .+ B*u
+measurement(x,u,p,t) = C*x
+
+# sss = SS(m,derivatives=false) |> collect
+# I(m.timings.nVars)[m.timings.past_not_future_and_mixed_idx,:] * [sss...,0,0]
+
+T    = 200 # Number of time steps
+kf   = KalmanFilter(A, B, C, 0, eye(nx), eye(ny), d0)
+ukf  = UnscentedKalmanFilter(dynamics, measurement, eye(nx), eye(ny), d0; ny, nu)
+x,u,y = LowLevelParticleFilters.simulate(kf,T,du) # Simuate trajectory using the model in the filter
+
+using SparseArrays, LinearOperators, BenchmarkTools
+n = 20
+p = 10
+x = sprand(p,n^3,.1)
+kk = sprand(n,n,.2) 
+kkLinOp = kk' |> LinearOperator
+
+@benchmark kron(kron(kron(kkLinOp, kkLinOp), kkLinOp), I(p)) * vec(x)
+@benchmark x * kron(kron(kk,kk),kk)
+
+sparse(reshape(kron(kron(kron(kk', kk'), kk'), I(p)) * vec(x),p,n^3))
+
+findnz(sparse(reshape(kron(I(p),kron(kron(kk,kk),kk)) * vec(x),p,n^3)))
+
+x * kron(kron(kk,kk),kk)
+findnz(x * kron(kron(kk,kk),kk))