Use Kalman Filtering to find the (Gaussian) loglikelihood, or for forecasting or smoothing.
1 2 3 4 5 6 7 8  KalmanLike(y, mod, nit = 0L, update = FALSE)
KalmanRun(y, mod, nit = 0L, update = FALSE)
KalmanSmooth(y, mod, nit = 0L)
KalmanForecast(n.ahead = 10L, mod, update = FALSE)
makeARIMA(phi, theta, Delta, kappa = 1e6,
SSinit = c("Gardner1980", "Rossignol2011"),
tol = .Machine$double.eps)

y 
a univariate time series. 
mod 
a list describing the statespace model: see ‘Details’. 
nit 
the time at which the initialization is computed.

update 
if 
n.ahead 
the number of steps ahead for which prediction is required. 
phi, theta 
numeric vectors of length ≥ 0 giving AR and MA parameters. 
Delta 
vector of differencing coefficients, so an ARMA model is
fitted to 
kappa 
the prior variance (as a multiple of the innovations variance) for the past observations in a differenced model. 
SSinit 
a string specifying the algorithm to compute the

tol 
tolerance eventually passed to 
These functions work with a general univariate statespace model with state vector a, transitions a < T a + R e, e ~ N(0, kappa Q) and observation equation y = Z'a + eta, eta ~ N(0, kappa h). The likelihood is a profile likelihood after estimation of kappa.
The model is specified as a list with at least components
T
the transition matrix
Z
the observation coefficients
h
the observation variance
V
RQR'
a
the current state estimate
P
the current estimate of the state uncertainty matrix Q
Pn
the estimate at time t1 of the state
uncertainty matrix Q (not updated by KalmanForecast
.
KalmanSmooth
is the workhorse function for tsSmooth
.
makeARIMA
constructs the statespace model for an ARIMA model,
see also arima
.
The statespace initialization has used Gardner et al's method
(SSinit = "Gardner1980"
), as only method for years. However,
that suffers sometimes from deficiencies when close to nonstationarity.
For this reason, it may be replaced as default in the future and only
kept for reproducibility reasons. Explicit specification of
SSinit
is therefore recommended, notably also in
arima()
.
The "Rossignol2011"
method has been proposed and partly
documented by Raphael Rossignol, Univ. Grenoble, on 20110920 (see
PR#14682, below), and later been ported to C by Matwey V. Kornilov.
It computes the covariance matrix of
(X_{t1},...,X_{tp},Z_t,...,Z_{tq})
by the method of difference equations (page 93 of Brockwell and Davis),
apparently suggested by a referee of Gardner et al (see p.314 of
their paper).
For KalmanLike
, a list with components Lik
(the
loglikelihood less some constants) and s2
, the estimate of
kappa.
For KalmanRun
, a list with components values
, a vector
of length 2 giving the output of KalmanLike
, resid
(the
residuals) and states
, the contemporaneous state estimates,
a matrix with one row for each observation time.
For KalmanSmooth
, a list with two components.
Component smooth
is a n
by p
matrix of state
estimates based on all the observations, with one row for each time.
Component var
is a n
by p
by p
array of
variance matrices.
For KalmanForecast
, a list with components pred
, the
predictions, and var
, the unscaled variances of the prediction
errors (to be multiplied by s2
).
For makeARIMA
, a model list including components for
its arguments.
These functions are designed to be called from other functions which check the validity of the arguments passed, so very little checking is done.
Durbin, J. and Koopman, S. J. (2001) Time Series Analysis by State Space Methods. Oxford University Press.
Gardner, G, Harvey, A. C. and Phillips, G. D. A. (1980) Algorithm AS154. An algorithm for exact maximum likelihood estimation of autoregressivemoving average models by means of Kalman filtering. Applied Statistics 29, 311–322.
R bug report PR#14682 (20112013) https://bugs.rproject.org/bugzilla3/show_bug.cgi?id=14682.
arima
, StructTS
. tsSmooth
.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15  ## an ARIMA fit
fit3 < arima(presidents, c(3, 0, 0))
predict(fit3, 12)
## reconstruct this
pr < KalmanForecast(12, fit3$model)
pr$pred + fit3$coef[4]
sqrt(pr$var * fit3$sigma2)
## and now do it year by year
mod < fit3$model
for(y in 1:3) {
pr < KalmanForecast(4, mod, TRUE)
print(list(pred = pr$pred + fit3$coef["intercept"],
se = sqrt(pr$var * fit3$sigma2)))
mod < attr(pr, "mod")
}

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.
All documentation is copyright its authors; we didn't write any of that.