eflik {KFAS} | R Documentation |
Function eflik computes log-likelihood function of univariate exponential family state-space model via simulation.
eflik(yt, Zt, Tt, Rt, Qt, a1, P1, P1inf, dist=c("Poisson", "Binomial", "Negative binomial"), offset=1,nsim=1000,tol=1e-10)
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. |
tol |
Tolerance parameter. Smallest covariance/variance value not counted for zero in diffuse phase. Default is 1e-10. |
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.
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. |
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
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 } #not run #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 } #not run #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