Continuous-Discrete Filtering using the Zakai Equation: Smooth - - PowerPoint PPT Presentation

Continuous-Discrete Filtering using the Zakai Equation: Smooth Likelihood Surface

Hermann Singer Department of Economics FernUniversität in Hagen Statistische Woche Trier September 2019

October 16, 2019

1/21

Continuous Time State Space Model

dY (t) = f (Y (t), t)dt + G(Y (t), t)dW (t) dZ(t) = h(Y (t), t)dt + dV (t) sampled measurements: zi = h(Y (ti), ti) + ǫi Goal: Optimal Filtering and Maximum Likelihood Estimation Wiener process W (t) Itô stochastic differential equations sampled measurements ˙ Z(ti) = zi, ρ(t)/dt = R(t)

2/21

(some) Solution Methods

Kalman Filtering (sequential)

Moment based methods

Taylor expansion: EKF, SNF, HNF Numerical integration: UKF, GHF, Smolyak sparse grid

PDE based methods: Stratonovich-Kushner and Duncan-Mortensen-Zakai (DMZ) equation Exact filters: Daum, Benes Particle Filters: Sequential Monte Carlo

Markov Chain Monte Carlo (nonsequential)

Simulated likelihood Bayesian approaches

3/21

Nonlinear bistable diffusion: Ginzburg-Landau model

dY = −[αY + βY 3]dt + σdW (t) = −∇Φ(Y ) + σdW (t) zi = Y (ti) + ǫi

• Figure: Simulated data (left) and extended Kalman filter (right).

4/21

Likelihood surface: Particle filter and Gauß-Hermite filter

• { }
• Likelihood surface particle filter {sample size} = {1000}

Figure: SIR particle filter (mean and SD) and trajectories (top, right), likelihood and score for β (bottom). Increment dβ = 0.0025.

5/21

Economic example: Equilibrium model: Herings (1996)

Potential Φ(y) ∼ α

2 y2 + β 4 y4

6/21

State estimation: Continuous-discrete Kalman Filter

time update: ti ≤ t < ti+1 (Fokker-Planck equation) ∂tp(y, t|Z i) = F(y, t)p(y, t|Z i) measurement update: t = ti+1 (Bayes formula) p(yi+1, ti+1|zi+1, Z i) = p(zi+1, ti+1|yi+1, Z i)p(yi+1, ti+1|Z i) p(zi+1, ti+1|Z i) filter density p(y, t|Z t) Fokker-Planck operator F(y, t) = −∂αfα + 1

2∂α∂βΩαβ

7/21

Recursive likelihood

p(zi+1|Z i; ψ) =

• p(zi+1|yi+1, Z i)p(yi+1|Z i)dyi+1

:=

• u(yi+1|zi+1, Z i)dyi+1

≈

• l

wl u(yi+1,l|zi+1, Z i) unnormalized filter density u(y, t|Z t) numerical integration using quadrature formulas measurements up to time ti: Z i = {Z(s) | s ≤ ti}

8/21

Continuous time filtering: DMZ equation

SPDE: Zakai (1969) ∂tu(y, t|Z t) = [F + h′ρ−1( ˙ Z − h/2)] ◦ u(y, t|Z t) = [F(y, t) + M(y, t)] ◦ u(y, t|Z t) measurement precision ρ−1(y, t) =

i π(t − ti)(Ridt)−1

measurement density p(dZ(t)|y, Z t) ∝ exp

• − 1

2(dZ − hdt)′(ρdt)−1(dZ − hdt)

• dZ ◦ u : symmetrized product

9/21

Stochastic Representation: Feynman-Kac Formula

u(y, t|Z t) = E

• e

t

t0 M(Y (τ),τ)dτδ(y − Y (t))

• Z t

use Lie -Trotter and Zassenhaus formula e(F+M)δtδ(y − y′) ≈ eMδteFδtδ(y − y′) = eMδtp(y, t + δt|y′, t) dY (t) = f (Y , t)dt + G(Y , t)dW (t), Y (t0) ∼ p(y, t0|Z t0) Dirac delta function limn→∞

• δn(x)φ(x)dx = φ(0)

10/21

Importance Sampling: Backward DMZ Equation

time reversal c(x, s) = u(x, T − s), s ≤ T ∂sc + Lc + (M + v)c = terminal condition c(x, T) = h(x) = u(x, 0) c(x, s) = E

• e

T

s (M+v)(X(τ),τ)dτh(X(T))

• X(s) = x, Z T−s

dX(τ) = ˜ f (X, T − τ)dt + G(X, T − τ)dW (τ), X(s) = x backward operator L = [−fα + (∂βΩαβ)]∂α + 1

2Ωαβ∂α∂β

scalar potential v = −(∂αfα) + 1

2(∂α∂βΩαβ)

11/21

Simulation of Backward DMZ Equation

Stochastic representation c(x, s) = E

• e

T

s (M+v)(X(τ),τ)dτh(X(T))

• X(s) = x
• dX(τ)

= ˜ f (X, T − τ)dt + G(X, T − τ)dW (τ) X(s) = x Importance sampling: drift correction (Milstein; 1995) Ω(X, T − τ)∇ log u(X, T − τ) approximate filter solution ˆ u(X, T − τ): EKF, GHF, UKF or particle filter

12/21

Ginzburg-Landau Model: Forward and backward simulation

• { _= *= { }={ } = _= = }
• Figure: Estimated filter density (top, right), backward simulation

(bottom) and forward simulation (top, left)

13/21

Ginzburg-Landau Model: Forward and backward simulation

• { _= *= { }={ } = _= = }
• Figure: Estimated filter density (top, right), backward simulation

(bottom) and forward simulation (top, left)

14/21

Likelihood surface: Particle filter and Zakai Equation

• { }
• Likelihood surface particle filter {sample size} = {1000}
• { }
• Likelihood surface Zakai filter {sample size,beta,alpha,dx} = {20, 1.5, 1.5, {0.25}}

Figure: Likelihood for SIR particle filter (top) and ZKF (Riemann), GHF, TKF (Riemann). Increment dβ = 0.0025.

15/21

Likelihood surface: Zakai Equation (UT)

• { }
• Likelihood surface Zakai filter {sample size,beta,alpha,method} = {20, 1.5, 1.5, {UT, 2, _}}
• { }
• Likelihood surface Zakai filter {sample size,beta,alpha,method} = {100, 1.5, 1.5, {UT, 2, _}}

Figure: Likelihood for ZKF (unscented transform UT), GHF, TKF. N = 20, 100. Increment dβ = 0.0025.

16/21

Conclusions

Continuous-discrete filtering with continuous time measurement equation Feynman-Kac representation

• f backward Zakai equation

Variance reduced simulation

• f unnormalized filter density

at supporting points No resampling required Smooth likelihood approximation using quadrature formulas at supporting points

17/21

References

Blankenship, G. and Baras, J. (1981). Accurate evaluation of stochastic Wiener integrals with applications to scattering in random media and to nonlinear filtering, SIAM Journal on Applied Mathematics 41(3): 518–552. Daum, F. and Huang, J. (2016). A plethora of open problems in particle flow research for nonlinear filters, Bayesian decisions, Bayesian learning, and transport, Signal Processing, Sensor/Information Fusion, and Target Recognition XXV, Vol. 9842, International Society for Optics and Photonics, p. 98420I. Herings, J. (1996). Static and Dynamic Aspects of General Disequilibrium Theory, Kluwer, Boston, London, Dordrecht. Hürzeler, M. and Künsch, H. R. (2001). Approximating and maximising the likelihood for a general state-space model, Sequential Monte Carlo methods in practice, Springer, pp. 159–175. Kantas, N., Doucet, A., Singh, S. S., Maciejowski, J., Chopin, N. et al. (2015). On Particle Methods for Parameter Estimation in State-Space Models, Statistical Science 30(3): 328–351.

17/21

Lemos, J. M., Costa, B. A. and Rocha, C. (2018). A Fokker-Planck approach to joint state-parameter estimation, IFAC-PapersOnLine 51(15): 389–394. Malik, S. and Pitt, M. K. (2011). Particle filters for continuous likelihood evaluation and maximisation, Journal of Econometrics 165(2): 190–209. Milstein, G. N. (1995). Numerical integration of stochastic differential equations, Vol. 313, Springer Science & Business Media (1988 in Russian). Mitter, S. K. (1982). Nonlinear Filtering of Diffusion Processes: A Guided Tour, Advances in Filtering and Optimal Stochastic Control, Springer, pp. 256–266. Singer, H. (2014). Importance Sampling for Kolmogorov Backward Equations, Advances in Statistical Analysis 98: 345–369. Zakai, M. (1969). On the optimal filtering of diffusion processes, Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete 11(3): 230–243.

18/21

Operator splitting

Lie –Trotter formula lim

n→∞[eAt/neBt/n]n

= e(A+B)t Zassenhaus formula eλ(A+B) = eλAeλBeλ2C2eλ3C3.... C2 =

1 2[B, A]

C3 =

1 3[C2, A + 2B]

e(A+B) ≈

• eA/neB/neC2/n2eC3/n3...eCm/nmn

18/21

Stratonovich calculus

dZ(t)u(y, t) = dZ(t) ◦ u(y, t) − 1

2h(y, t)u(y, t)dt

DMZ equation in Itô-form du(y, t|Z t) = [F(y, t)dt + h′(y, t)ρ−1(t)dZ(t)]u(y, t|Z t) symmetrized product dZ(t) ◦ u(y, t) := dZ(t)¯ u(y, t) ¯ u(y, t) :=

1 2[u(y, t) + u(y, t + dt)]

u(y, t) = ¯ u(y, t) − 1

2du(y, t)

19/21

Potential Φ(y) = α

2y 2 + β 4y 4,

drift f (y) = −∇Φ

• 10
• 5

5 10

• 20
• 10

10 20

• 10
• 5

5 10 0.05 0.1 0.15 0.2 0.25 0.3

Figure: Left: Potential as a function of y for parameter values α = −3, −2, ..., 1. Right: Stationary density pstat ∝ exp[−(2/σ2)Φ(y)].

20/21

Importance sampling: Kolmogorov Backward Equation

∂sc(x, s) + L(x, s)c(x, s) + v(x, s)c(x, s) = terminal condition c(x, T) = h(x) solution c(x, s) = E

• e

T

s v(Y (τ),τ)dτh(X(T))

• X(s) = x
• dX(t) = f (X, t)dt + G(X, t)dW (t), X(s) = x

importance sampling: drift correction Ω(x, s)∇ log c(x, s) (Milstein; 1995) backward operator L = fα∂α + 1

2Ωαβ∂α∂β

21/21