### Kalman Smoother and EM for Stochastic Constant Growth Model

#### by allenlu2007

以下的例子為 diffusion approximation for a stochastic exponential growth model. 可以用 linear space model

B 是平均 growth rate. A = C = 1.

Initial condition t=1

再來結合 Markov 的 property (1)

and observed output only depends on the corresponding states (2)

We can get a very simple joint pdf of all Xi and Yi in Kalman filter (or Hidden Markov Model) as!!

以下是有趣的部份:

Where only Yi are given.

記住 Xi 是 hidden variables. 因些我們可以(需要)用 EM algorithm 來 optimize

## EM Algorithm

由 Shumway and Stoffer (1982) 所建議，基本上有四個步驟。

(0) 計算 initial parameter estimates

(1) Kalman smoother: generate an estimate of xt by estimating E(xt | Y1, … YT) given above initial values.

(2) Define Q function as below. Update R, Q, B, Pi, V by maximizing the Q function.

(3) Check if si function is converged and repeat (1)

以下為分解動作

### Step 0: Initial Parameter Estimates

By Holmes and Fagan (2002)

### Step 1: Kalman Smoother (Kalman-Rauch Recursion)

(a) 先用 Kalman filter forward pass to estimate Gaussian pdf of p(xt | y1, .. yt).

Gaussian pdf 由 mean and variance 完全決定。Mean (xt = E(xt | y1, .. yt) 和 variance (Vt) 結果如下

Vt 在 prediction phase 時會增加 (by Q), 而在 update phase 時會減小 (by Kt Vt^(t-1))

(b) 再用 Rauch backward pass to estimate p(xt| y1, .. yT). 同樣 mean (xt) and variance (Vt) 結果如下

最後一項的右手第二項一般是負的，所以整體的 variance 可以再減少。

(c) 再一個 backward pass 是計算

### Step 2: EM Algorithm

接下來定義 E-step and M-step. E-step 如下

利用在開始所推導的 log P 公式，其實我們不必算出 psi function.

重點在於 taking partial derivative of psi w.r.t. each parameters, set the derivative to zero and solve the maximizing parameter:

(好像少了一個 1/T log |R| 的導數? and what is C?)

### Step 3: 檢查收歛性

Matlab code example

clear all

close all

A = 1

C = 1

B = 0.2

D = 0

t = [0:200];

Q = 0.02;

R = 1;

x0 = 15;

sys = ss(A, B, C, D, -1);

u = ones(length(t), 1); % external force B*u

v = sqrt(Q)*randn(length(t), 1); %state noise

w = sqrt(R)*randn(length(t), 1); %measurement noise

y = lsim(sys, u+v, t, x0);

logdata = y + w;

plot(logdata)

[BB, QQ, RR, initx, V1, LL] = PVAKalman(logdata’)

BB = 0.2

RR = 0.9748

QQ = 0.0013

initx = 15.1161

## EM Algorithm 2:

考慮另一個 model

A, B, and C 是已知，但是 Q and R 未知: