EMalgo: Estimation of variance matrices in a Gaussian state space...
In ClausDethlefsen/sspir: State Space Models in R

Description Usage Arguments Details Value Author(s) See Also Examples

Estimates variance matrices of the observation- and latent process in a Gaussian state space model given as input, using the EM algorithm.

EMalgo(	ss,
        epsilon   = 1e-06,
        maxiter   = 50,
        Vstruc    = function(V) { return(V) },
        Wstruc    = function(W) { return(W) }
      )

`ss`	an object of class `SS`.
`epsilon`	a (small) positive numeric giving the tolerance of the maximum relative differences of `Vmat` and `Wmat` between iterations.
`maxiter`	a positive integer giving the maximum number of iterations to run.
`Vstruc`	a function specifying the structure of the variance matrix of the observation model if such is available.
`Wstruc`	a function specifying the structure of the variance matrix of the state model if such is available.

The EM algorithm is a two-step iterative maximisation algorithm used to estimate unknown and constant covariance matrices V and W in a Gaussian state space model. Let the state space model be given as in SS with constant covariance matrices.

The covariance matrix V is estimated by maximisation of log(p(Y_t|θ_t, V)) conditional on all information, denoted D_n. This is equivalent to minimise the expectation of

log(V) + n^-1 Σ (Y_t - μ_t)^2V^-1

conditional on D_n.

This yields

E[log(V) + n^-1 Σ (Y_t - μ_t)^2 V^-1 | D_n]

log(V) + n^-1 Σ [ trace( V^-1 F_t^T s.C F_t ) + (Y_t - F_t^T s.m_t)^2]

log(V) + trace[ (nV)^-1 Σ ( F_t^T s.C F_t + (Y_t - F_t^T s.m_t)^2)].

This is minimised for

V^* = n^-1 Σ [ F_t^T s.C F_t + (Y_t - F_t^T s.m_t)^2].

Similarly W is estimated by minimising the expectation of

log|W| + n^-1 Σ (θ_t - G_tθ_t-1)(θ_t - G_tθ_t-1)^T W^-1

conditional on D_n.

This yields

E[log|W| + n^-1 Σ (θ_t - G_tθ_t-1)(θ_t - G_tθ_t-1)^T W^-1 | D_n]

log|W| + n^-1 Σ [ trace( W^-1 L_t ) + (s.m_t - G_t s.m_t-1)(s.m_t - G_t s.m_t-1)^T]

log|W| + trace[ (nW)^-1 Σ ( L_t + (s.m_t - G_t s.m_t-1)(s.m_t - G_t s.m_t-1)^T)].

This is minimised for

W^* = n^-1 Σ [ L_t + (s.m_t - G_t s.m_t-1)(s.m_t - G_t s.m_t-1)^T],

where s.m_t, s.C_t denotes smoothed means and covariances and L_t=Var[θ_{t}-G_tθ_{t-1}|D_n].

`ss`	the output from `smoother`.
`Vmat.est`	the estimate of the observation variance matrix, which is provided if the input of `Vstruc` is of class `function`, otherwise as input of `Vmat`.
`Wmat.est`	the estimate of the observation variance matrix, which is provided if the input of `Wstruc` is of class `function`, otherwise as input of `Wmat`.
`loglik`	maximum value of log likelihood function for each iteration.
`iteration`	number of iterations upon convergence.

Anette Luther Christensen

kfilter, smoother, recursion.

## Formulates Gaussian state space model:
## Trend: local linear
## Seasonal variation: sine function with frequency = 1
m1 <- SS( Fmat = function(tt, x, phi) {
            Fmat      <- matrix(NA, nrow=3, ncol=1)
            Fmat[1,1] <- 1
            Fmat[2,1] <- cos(2*pi*tt/12)
            Fmat[3,1] <- sin(2*pi*tt/12)
            return(Fmat)
          },
          Gmat = function(tt, x, phi) {
            matrix(c(1,0,0,0,1,0,0,0,1), nrow=3)
          },
          Wmat = matrix(c(0.01,0,0,0,0.1,0,0,0,0.1), nrow=3),
          Vmat = matrix(1),
          m0   = matrix(c(0,0,0), nrow=1),
          C0   = matrix(c(1,0,0,0,0.001,0,0,0,0.001), nrow=3, ncol=3)
        )

## Simulates 100 observation from m1
m1 <- recursion(m1, 100)

## Specifies the correlation structure of W:
Wstruc <- function(W){	
            W[1,2:3] <- 0
            W[2:3,1] <- 0
            W[2,2]   <- (W[2,2]+W[3,3])/2
            W[3,3]   <- W[2,2]
            W[2,3]   <- 0
            W[3,2]   <- 0
            return(W)
          }

## Estimstes variances and covariances of W by use of the EM algorithm
estimates <- EMalgo(m1, Wstruc=Wstruc)
estimates$ss$Wmat

## Plots estimated model
plot(estimates$ss)