Spectral LogLikelihood Function and Derivatives
Description
These functions evaluate the negative of the spectral loglikelihood function of a linear Gaussian state space model and its first and second order derivatives.
Usage
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19  mloglik.fd(x, model,
barrier = list(type = c("1", "2"), mu = 0), inf = 99999, xreg)
mcloglik.fd(x, model, xreg = NULL,
barrier = list(type = c("1", "2"), mu = 0), inf = 99999)
mloglik.fd.deriv(model, xreg = NULL,
gradient = TRUE, hessian = TRUE, infomat = TRUE, modcovgrad = TRUE,
barrier = list(type = c("1", "2"), mu = 0),
version = c("2", "1"))
mcloglik.fd.deriv(model, xreg = NULL,
gradient = TRUE, hessian = TRUE, infomat = TRUE)
mloglik.fd.grad(x, model, xreg = NULL,
barrier = list(type = c("1", "2"), mu = 0),
inf)
mcloglik.fd.grad(x, model, xreg = NULL, inf, barrier)

Arguments
x 
a numeric vector containing the parameters of the model.
This is an auxiliary argument so that this function can be used as input
to 
model 
an object of class 
xreg 
optional list containing constant terms. See details. 
barrier 
a list defining a barrier term to penalize parameter values close to the bounds

inf 
a numeric indicating the value to be returned if the value of the loglikelihood
function happens to be 
gradient 
logical. If 
hessian 
logical. If 
infomat 
logical. If 
modcovgrad 
logical. If 
version 
a character indicating whether implementation 
Details
The spectral loglikelihood of a linear Gaussian state space model is given by (Harvery, 1989 Section 4.3):
logLik = 0.5 log(2π)  0.5 ∑_{j=0}^{n1} log\,g(λ[j])  π ∑_{j=0}^{n1} I(λ[j])/g(λ[j])
where λ[j] is a frequency defined as λ[j] = 2π j/n; I(λ[j]) is the periodogram at frequency λ[j] and g(λ[j]) is the spectral generating function of the model at frequency λ[j].
The derivation of the spectral likelihood function defined above relies on the assumption that the process is circular (its covariance matrix is circulant). If the process is not circular the value of the likelihood is an approximation.
First and second order derivatives are computed by means of their analytical expressions. The first order derivatives of the spectral loglikelihood with respect to parameter θ are given by:
(d\,logLik)/(d\,θ) = 0.5 ∑_{j=0}^{n1} ( (2π I(λ[j])) / (g(λ[j]))  1 ) (1/g(λ[j])) d g(λ[j])/d θ
Second order derivatives are given by:
(d^2\, logLik)(d\,θ θ') = ∑_{j=0}^{n1} ( (2π I(λ[j]))/(g(λ[j]))  1 ) 1/(2g(λ[j])) (d^2 g(λ[j]))/(d θ d θ') 
2 ∑_{j=0}^{n1} ( (4π I(λ[j]))(g(λ[j]))  1 ) ( \frac{1}{2g(λ[j])} )^2 (d g(λ[j]))/(d θ) (d g(λ[j]))/(d θ')
The argument x
is an auxiliary vector that is necessary in some contexts.
For example, the input to function optim
must contain
as first argument the vector of parameters where optimization is performed.
If it is not required or is redundant information contained in model@pars
it can be set to NULL
.
The functions mcloglik.fd
, mcloglik.fd.deriv
and
mcloglik.fd.grad
use the expressions for the spectral loglikelihood
function where the parameter specified in model@cpar
is concentrated
out of the likelihood function.
For further information about the barrier term see
Bounds on parameters and barrier term
in the details section in maxlik.fd.scoring
.
Arguments inf
and barrier
are not used by mloglik.fd.grad
and mcloglik.fd.grad
but they are needed in maxlik.fd.optim
,
where this function is passed as the gradient to be used by optim
including the arguments inf
and barrier
.
Argument xreg
.
It is an optional list of constant terms. It is used by
maxlik.fd.optim
when analytical derivatives are employed and
by maxlik.fd.scoring
. It avoids computing some constant terms each
time the function mloglik.fd.grad
is called.
The list xreg
should contain
an element called dxreg
, the external regressors differeneced by means
of the differencing filter that renders stationarity in the model and
the element fft.xreg
, the Fourier transform of each regressor
in dxreg
.
The list xreg
is not used by mloglik.fd
. It is necessary
to define this argument in the prototype of the function because when this function
is passed to optim
along with
mloglik.fd.grad
, the argument xreg
is passed to mloglik.fd
when it is defined in optim
as an argument to be passed
to mloglik.fd.grad
.
Argument xreg
is not currently implemented in functions with
concentration of a parameter, mloglik.fd
, mcloglik.fd.deriv
and mcloglik.fd.grad
.
Note: modcovgrad
is not available when external regressors are defined
in the input model, model
.
Value
mloglik.fd
returns a numeric value of the negative
of the spectral loglikelihood evaluated at the parameter values defined in
the input model model
or at x
if this argument is not NULL
.
If the value happens to be NA
or nonfinite the value of argument inf
is returned. This function is suited to be passed as the objective function to
optim
.
mloglik.fd.deriv
returns a list containing a vector of the first order derivatives
of the negative of the spectral likelihood and a matrix for the second order
derivatives. Those derivative terms that are not requested by setting the
corresponding argument to FALSE
are set to NULL
in the output list.
mloglik.fd.grad
returns a numeric vector containing the gradient.
This function is suited to be passed as the gradient function to
optim
.
mcloglik.fd
, mcloglik.fd.deriv
and mcloglik.fd.grad
return the value of the same information as the other functions but for the
concentrated likelihood function.
Note
mcloglik.fd.deriv
is not currently implemented for model with nonnull
model@transPars
or with a barrier term.
References
Harvey, A. C. (1989). Forecasting, Structural Time Series Models and the Kalman Filter. Cambridge University Press.
See Also
barrier.eval
,
logLik
,
maxlik.fd
,
stsm
.
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26  # define the local level model for Nile time series
pars < c("var1" = 11000, "var2" = 1700)
m < stsm.model(model = "locallevel", y = Nile, pars = pars)
# 'mloglik.fd' returns the negative of the loglikelihood function
mloglik.fd(model = m)
# 'logLik' returns the value of the loglikelihood function
logLik(object = m, domain = "frequency")
# compare analytical and numerical derivatives
# more tests in file 'testderivativesmloglikfd.R' in the
# folder 'inst' of the source package
system.time(da < mloglik.fd.deriv(m, gradient = TRUE, hessian = TRUE))
dgn < numDeriv::grad(func = mloglik.fd, x = m@pars, model = m)
dhn < numDeriv::hessian(func = mloglik.fd, x = m@pars, model = m)
all.equal(as.vector(da$gradient), dgn)
all.equal(da$hessian, dhn)
# the same as above for the local level plus seasonal model and
# a sample simulated series
data("llmseas")
m < stsm.model(model = "llm+seas", y = llmseas)
system.time(a < mloglik.fd.deriv(model = m, gradient = TRUE, hessian = TRUE))
system.time(g < numDeriv::grad(func = mloglik.fd, x = m@pars, model = m))
system.time(h < numDeriv::hessian(func = mloglik.fd, x = m@pars, model = m))
all.equal(a$gradient, g, check.attributes = FALSE)
all.equal(a$hessian, h, check.attributes = FALSE)
