KFAS (version 0.6.1)

eflik: Log-likelihood of exponential family state-space model.

Description

Function eflik computes log-likelihood function of univariate exponential family state-space model via simulation.

Usage

eflik(yt,  Zt,  Tt,  Rt,  Qt,  a1,  P1,  P1inf,  dist=c("Poisson",
"Binomial",  "Negative binomial"),  offset=1,  nsim=1000, maxiter=50)

Arguments

yt

Matrix, array or vector of observations. Note that yt is univariate.

Zt

System matrix or array of observation equation.

Tt

System matrix or array of transition equation.

Rt

System matrix or array of transition equation.

Qt

Variance matrix or array of disturbance terms eta_t.

a1

Initial state vector.

P1

Variance matrix of a1. In diffuse case P1star, the non-diffuse part of P1.

P1inf

Diffuse part of P1.

dist

Distribution of yt.

offset

Vector of length n. See details.

nsim

Number of simulations.

maxiter

Maximum number of iterations used in linearisation.

Value

• List with output from Kalman smoother, when model is approximated with gaussian distribution g(y|theta). Note that Ht is H*_t. List also contains following items:
• ytilde
  y* of approximating gaussian model.
• theta
  Z_t * alpha_t of approximating gaussian model.
• likp
  Value of logL_p.
• offset
• dist
  Distribution of yt.
• lik
  Value of log-likelihood function.

Details

Function computes log-likelihood of exponential family state-space model. It is very recommended to use estimates gained from function lik0 as initial values.

logL_p(psi) = logL_g(psi) + logE*_g(w|y),

where logL_g(psi) is log-likelihood of approximate gaussian distribution and logE*_g(w|y) is a Monte-Carlo approximation of E_g[p(y|theta)/g(y|theta) | y]. For details, see J. Durbin and S.J. Koopman (1997).

The general state space model for exponential family is given by

p(y_t|theta_t) = exp[theta'_t * y_t - b_t(theta_t) + f_t(y_t)] (observation equation)

alpha_t+1 = T_t * alpha_t + R_t * eta_t (transition equation)

where theta_t = Z_t * alpha_t and eta_t ~ N(0,Q_t).

Approximating gaussian model is given by

y*_t = Z_t * alpha_t + eps_t (observation equation)

alpha_t+1 = T_t * alpha_t + R_t * eta_t (transition equation)

where eps_t ~ N(0,H*_t) and eta_t ~ N(0,Q_t).

If yt is Poisson distributed, parameter of Poisson distribution is offset*lambda and theta = log(lambda).

If yt is from binomial distribution, offset is a vector specifying number of trials at times 1,...,n, and theta = log[pi_t/(1-pi_t)], where pi_t is the probability of success at time t.

In case of negative binomial distribution, offset is vector of specified number of successes wanted at times 1,...,n, and theta = log(1-pi_t).

Note that this function works only for univariate observations.

References

Durbin J. and Koopman, S.J. (1997). Monte Carlo Maximum Likelihood Estimation for Non-Gaussian State Space Models, Biometrica, Vol. 84, No. 3. Koopman, S.J. and Durbin J. (2001). Time Series Analysis by State Space Methods. Oxford: Oxford University Press

Examples

kk <- c(396, 399, 403, 434, 340, 519, 558, 566, 591, 566, 574, 646, 644,
646, 665, 693, 773, 834, 910, 1035, 1002, 1161, 1056, 1097, 1094, 1042,
1194, 1316, 1246, 1503, 1428, 1477, 1490, 1465, 1560, 1860, 2008, 2020, 2167)

vlk <- c(4623785, 4606307, 4612124, 4639657, 4666081, 4690574, 4711440,
4725664, 4738902, 4752528, 4764690, 4779535, 4799964, 4826933, 4855787,
4881803, 4902206, 4918154, 4932123, 4946481, 4964371, 4986431, 5013740,
5041992, 5066447, 5088333, 5107790, 5124573, 5139835, 5153498, 5165474,
5176209, 5188008, 5200598, 5213014, 5228172, 5246096, 5266268, 5288720)


n<-39
Zt<-array(c(1,0),c(1,2))
Tt<-diag(c(1,1))
Tt[1,2]<-1
Rt <- diag(2)
a1 <- array(0,dim=2)
P1 <- diag(0,2)
P1inf <-diag(1,2)

likfn<-function(pars,yt,Zt,Tt,Rt,a1,P1,P1inf,dist,offset=1)
{
Qt<-diag(exp(pars))
lik<-eflik0(yt,Zt,Tt,Rt,Qt,a1,P1,P1inf,dist,offset)$lik0
lik
}

opt <- optim(par=c(1,1), yt=array(kk,dim=c(1,39)), Zt=Zt, Tt=Tt, Rt=Rt,
a1=a1, P1=P1, P1inf=P1inf, offset=vlk, dist="Poisson", fn=likfn,
method="BFGS", control=list(fnscale=-1,trace=1,REPORT=1))

Qt<-diag(exp(opt$par))
out<-eflik0(yt=kk,Zt,Tt,Rt,Qt,a1,P1,P1inf,offset=vlk)
out$Qt*10000


likfn2<-function(pars,yt,Zt,Tt,Rt,a1,P1,P1inf,dist,offset=1)
{
Qt<-diag(exp(pars))
lik<-eflik(yt,Zt,Tt,Rt,Qt,a1,P1,P1inf,dist,offset,nsim=1000)$likp
lik
}

opt2 <- optim(par=opt$par, yt=array(kk,dim=c(1, 39)), Zt=Zt, Tt=Tt,
Rt=Rt, a1=a1, P1=P1, P1inf=P1inf, offset=vlk, dist="Poisson", fn=likfn2,
method="BFGS", control=list(fnscale=-1, trace=1, REPORT=1))

Qt<-diag(exp(opt2$par))
out<-eflik(yt=kk, Zt, Tt, Rt, Qt, a1, P1, P1inf, offset=vlk)
out$Qt*10000